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1. Introduction and Nomenclature 

The linear-quadratic-gaussian (LQG) compensator [1-3] has been developed to facilitate the 
design of control laws for multi-input multi-output (MIMO) systems. An LQG compensator mini- 
mizes a quadratic performance index and (under mild conditions) is guaranteed to yield an internally 
stable closed-loop system. Unfortunately, however, the minimal dimension of an LQG compensator 
is almost always equal to the dimension of the plant and can thus often violate practical implemen- 
tation constraints on controller order. This deficiency is especially highlighted when considering 
control-design for high-order systems such as flexible space structures. Hence, a very relevant area 
of research is the development of methodologies that will enable the design of optimal controllers 
whose dimension is less than that of the design plant (i.e., reduced-order controllers). 

Two main approaches have been developed to tackle the reduced-order design problem. The 
first approach attempts to develop approximations to the optimal reduced-order controller by re- 
ducing the dimension of an LQG controller [4-11]. These methods are attractive because they 
require relatively little computation and should be used if possible. Unfortunately, they tend to 
yield controllers that either destabilize the system or have poor performance as the requested con- 
troller dimension is decreased and/or the requested authority level is increased. Hence, if used in 
isolation, these methods do not yield a reliable methodology for reduced-order design. 

The second approach attempts to directly synthesize an optimal, reduced-order controller by 
a numerical optimization scheme [12-25]. Almost all of these schemes are parameter optimization 
approaches; that is, they represent the controller by some parameter vector and attempt to find 
the vector that optimizes the cost functional. Unfortunately, most of these schemes have only 
local convergence properties and hence have the potential of failure if the initial controller is not 
“close” to the optimal controller. One exception is the homotopy algorithm described in [20,25]. A 
homotopy allows an initial controller to be deformed gradually into the desired optimal controller 
by following a homotopy path. These schemes are particularly useful because they have global 
convergence properties . Hence, this algorithm does not require the initial controller to be near 
the optimal controller. The algorithm is based on solving a set of “optimal projection” equations 
[26,27] that are a characterization of the necessary conditions for optimal reduced-order control. 
Unfortunately, the algorithm has sublinear convergence properties and the convergence slows at 
higher authority levels and may fail. 

This volume describes a new homotopy algorithm for discrete-time systems which has been 
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implemented in MATLAB. The homotopy algorithm is based on a parameter optimization for- 
mulation. This algorithm shares the global convergence properties of the homotopy algorithm of 
[20,25] but potentially has quadratic or superlinear convergence rates. The results reported here 
may offer the foundation for a reliable approach to optimal, reduced-order controller design. 


Nomenclature 

Y>Z 
Y > Z 
z ij)Zi,j or 

It 

z* 

Z* 


tr Z 

mu 

vec(-) 


era 


p(»J) 

^mXn 

r(0 

On 

Z 


Y * Z 
Y/Z 


x = A\b 



Y — Z is nonnegative definite 

Y — Z is positive definite 
(t,j) element of matrix Z 
r x r identity matrix 

the group generalized inverse of the square matrix Z 
satisfying rank Z = rank Z 2 [28,29] 

the (unique) nonnegative definite square root of Z (Z^Z* = Z), 
Z = Z T > 0 

trace of square matrix Z 

absolute norm of matrix Z (||Z||yi = max.j | Zij |) 
the invertible linear operator defined such that 


vec(s) = [$7 s J"' 5 J] T > § € nt*’* 9 
where Sj € IR P denotes the j ih column of S. 


the m-dimensional column vector whose i th element 
equals one and whose additional elements are zeros, 
the m x n matrix whose (i,j) element equals one 

T 

and whose additional elements are zero (= eme^ )* 


the m x m matrix whose » th row has all unity 
elements and whose additional rows are zero. 

for the square matrix Z, Z is the identically 
dimensioned matrix defined by Zij = Zii. 

Hadamard product of Y and Z([yijZij]) 

(Y and Z must have identical dimensions.) 

matrix whose (i,j) element is yij/zij 
(Y and Z must have identical dimensions.) 

(MATLAB notation) 

x is the least squares solution to Ax = b 

m x m matrix having unity elements (i.e., N m< ij = 1) 

matrix whose t‘ h row is given by the row vector 


1-2 


October 1993 


GASD-HADOC 


Harris Corporation 


i T and whose additional rows are zero. (The size of 
the matrix is understood from the context) 

matrix whose j ih column is given by the column vector x and whose 
additional columns are zero. (The size of 
the matrix is understood from the context.) 

A: th row of the matrix Z 
(MATLAB notation) 

k ih column of the matrix Z 
(MATLAB notation) 

fc th row of the matrix XY Z 

fc th column of the matrix XY Z 

partitioned symmetric matrix whose (1,1), (2,2) and (1,2) matrix 
partitions are as given. 
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2. Optimal Reduced-Order Discrete-Time Dynamic Compensation 
Consider the discrete-time system 

x(k + 1) = Ax(k) + Bu(k) + w\(k) (2.1a) 

y(k) = Cx(k) + Du(k) + w 2 (k) (2.1 b) 

where x G IR"*, u G nt"*, y G IR' 1 '', w\ G IR n * is a white noise disturbance with covariance 
Vi > 0, w 2 G IR n * is white observation noise with covariance V 2 > 0, and Wj and w 2 have cross 
covariance V\ 2 G IR n * Xn *. If D = 0, we desire to design a fixed-order dynamic compensator, 


x c (k + 1) = A c x c (k) + B c y(k) (2-2 a) 

u(k) = -C c x c (k) - D c y(k), (2-2 b) 

or if D ^ 0, we desire to design a fixed-order dynamic compensator 

x c (k + 1) = A c x c (k ) + B c y(k) (2.3a) 

u(k) = —C c x c (k) (2.3b) 

which minimizes the steady-state performance criterion 

J(A c ,B c ,C c ,D c ) = lim E[x T (k)R 1 x(k) + 2x T (k)R l2 u(k) + u' T (k)R 2 u(k)] (2.4) 

k—+ OO 


where x c G IR n<: , n c < n x , R\ = Rj > 0, and R 2 = Rj > 0. We will call this problem 
the optimal reduced-order dynamic compensation problem for discrete-time systems. 

The closed-loop system corresponding to (2.1) and (2.2) or (2.1) and (2.3) can be expressed as 

x(k + 1) = -Az(A:) + w(k) (2-5) 


where 


*(*) = 


*(*) 
*c(*)J ’ 


w(k) = 


w\(k) — BD c w 2 (k) 

B c w 2 (k) 


- \ A - BD c C -BC c 

4 “ B c C A c - B c DC c 


( 2 . 6 ) 

(2.7) 


and it understood that either D of D c is identically zero in (2.61 and (2.7'). In addition, the cost 


(2.4) can be expressed as 


J(A c ,B c ,C c ,D c ) = hm E[x T (k)Rx(k)} + E[wJ(k)DjR 2 D c w 2 (k)] (2.8) 

k—* 00 
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where 


MI" f ,! l 

1&12 ^22 J 

(2.9) 

R n = R!- C T DjRj 2 - R 12 D c C + C T DjR 2 D c C 

(2.10a) 

R\ 2 = — R\ 2 C c + Dj R 2 C c 

(2.106) 

r 22 — CjR 2 C c . 

(2.10c) 


To guarantee that the cost J is finite and independent of initial conditions we restrict our 
attention to the set of stabilizing compensators, 


S c = {( A c ,B c ,C c ,D c ) : A is asymptotically stable}. 


( 2 . 11 ) 


Assume (A c , B c ,C c , D c ) € S c and define Q and P to be the closed-loop steady state covariance and 
its dual, i.e., 


Q = AQA t + V (2.12) 

P = A t PA + R ( 2 . 13 ) 


where 


and 



Viz' 

V 22 


( 2 . 14 ) 


V n =Vi- BD C V £ - V 12 DjB' 1 + BD c V 2 DjB T ( 2 . 15 a ) 

V n = V 12 Bj - BD c V 2 Bj ( 2 . 156 ) 

V 22 = B c V 2 Bj. ( 2 . 15 c ) 


Then, the cost can be expressed as 

J(A c ,B c ,C c ,D c ) = ti(QR) + ti(DjR 2 D c V 2 ). ( 2 . 16 ) 

Also, P and Q can be expressed in the partitioned forms 


P = 

\Pi i 

P T 

■*12 

Pn 

P 22 

, P„6lR niXni 

( 2 . 17 ) 

Q = 

Q ii 

m 

Qn 

, Qn € IR" 1 *"*. 

( 2 . 18 ) 
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Notice that Qu is the covariance of the plant states, Q 22 is the covariance of the compensator 
states and Q 22 is the cross-covariance of the plant and controller states. 

Since the value of J is independent of the internal realization of the compensator, in what follows 
we will further restrict our attention to minimal compensators. Hence, we define the admissable set , 

5+ = {(A c ,B c ,C c ,D c ) € S : (A C ,B C ) is controllable, (j4 c ,C c ) is observable}. (2.19) 

Note that S+ is an open set. 

Optimal Projection theory can be used to characterize all admissable extremals of the optimal 
reduced-order dynamic compensation problem for discrete-time systems. Before presenting the 
main theorems we present an important Lemma and some definitions which are useful in stating 
the main results of optimal projection theory. The lemma also gives many properties of the optimal 
projection solution (see Theorem 2.1). 

Lemma 2.1 [1]. Suppose Q E nt n * Xn ' and P E IR n * Xn * are symmetric and nonnegative 
definite and rank QP = n c . Then, the following statements hold: 

(i) QP is diagonalizable and has nonnegative eigenvalues. 

(ii) The n x x n x matrix 

r = QP{QP)* (2.20) 

is idempotent, i.e., r 2 = r (r is an oblique projection) and 

rank r = n c . (2-21) 

Thus, if t is given by (2.18), then there exists a nonsingular matrix W £ IR n * xn * such that 

T = W 0 W ~ 1 ' ( 2 ‘ 22 ^ 

(iii) There exists G,T £ Dl n<;Xnx and nonsingular M £ IR ncXn ‘ : such that 

QP = G T Mr (2.23) 

rG T = J nc . (2.24) 

(iv) If G,r and M satisfy property (iii) then 

rank G = rank f = rank M — n c (2.25) 
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(2.26) 

(2.27) 

(2.28) 


(QP)* = G T M~ 1 r 
t = G J r 

tG t = G T , Tr = r. 

(v) The matrices G, f and M satisfying property (iii) are unique except for a change of basis in 

IR"* , i.e., if G\ f and M' also satisfy property (iv), then there exists nonsingular T c G IR" 1 xnc 

such that G' = TjG, F = T~ l T, M' = T~ l MT c . Furthermore, all such M are diagonalizable 
with positive eigenvalues. 

(vi) Finally, if rank Q = rank P = rank QP = n c , there exists a nonsingular transformation 
W G nt n « xn, such that 

P = W~ T q ° W~ l (2.29) 

Q = W q J W T (2.30) 

where fl G IR n ‘ xnc is diagonal and nonsingular. In addition, 

P = r T P = Pt = t t Pt (2.3 1 ) 

Q = tQ = Qt t = tQt t . (2.32) 

Remark 2.1. The transformation W in statement (vii) meets the requirements of statement 

(ii). 

Definition 2.1. A triple (G,M,T) satisfying property (iii) of Lemma 2.1 is a projective 
factorization of QP. 

To optimize (2.8) subject to the constraint (2.12) we form the Lagrangian 

£{A c ,B c ,Cc, D C ,P,Q) = tr [QR + P(AQA T + V - Q) + DjR 2 D c V 2 } (2.33) 

where P G IR) n * +n ‘=) x ( n * +n ‘) is the Lagrange multiplier. The stationary conditions are then given 


by 



dC n 
— = 0, 
dp 

= o 
SQ 


(2.34) 

o 

II 

^ CO 

dL n 
dB c ~ ’ 

II 

© 

£ ^ 
II 

p 

(2.35) 
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Definition 2.2. A compensator ( A c ,B c ,C c ,D c ) is an extremal of the optimal reduced-order 
dynamic compensation problem for discrete-time systems if it satisfies the stationary conditions 
(2.32). 

Definition 2.3. A compensator ( A c ,B c ,C c ,D c ) is an admissible extremal of the optimal 
reduced-order dynamic compensation problem for discrete-time systems if it is an extremal and is 
also in S+ . 

Remark 2.2. The optimal (admissible) reduced-order dynamic compensator for discrete-time 
systems (if it exists) can be found by computing all admissible extremals. 

We can now state in the form of two theorems the basic result of Optimal Projection theory, 
namely a set of necessary conditions which characterize admissible extremals of the optimal fixed- 
order dynamic compensation problem. For convenience define 

P a = B t PA + Rj 2 , Q a = AQC T + V u (2. 36a, 4 ) 

Ri, a = R2+ B t PB, V 2 , a = V 2 + CQC t . (2.37a, 4 ) 

Theorem 2.1 [3]. Suppose D = 0 and (A c ,B c ,C c ,D c ) is an admissible extremal of the 
optimal reduced-order dynamic compensation problem for discrete-time systems. Then, there exist 
nonnegative-definite matrices P,Q,P and Q such that A c ,B c ,C c and D c are given by 


Ac = r(A - BR 2 ~iP a - Q a V 2 ~l + Q a V 2 }DR 2 \P a - BD c C)G T (2.38) 

B c = r(Q a V 2 -J + BD C ) (2.39) 

C c = (R 2 \P a + D c C)G t (2-40) 

D c = R 2 ^(B t PAQC t + Rj 2 QC T + B t PV 12 V 2 ~J (2.41) 


for some projective factorization ( G,M,T ) of Q,P such that the following conditions are satisfied: 
P =A t PA + R x - Pj R 2 \P a 

+ rl[{A - Q a V 2 ~lC) T P{A - QaV 2 -JC)+(P a + R2, a D c C) T R 2 ^(P a + R2,aD c C)r x (2.42) 
Q =AQA t + Vx - Q«V 2 ~lQl 

+ r ± [(A -BR 2 \P a )Q{A - BR^ a P a ) T + (Q a + BD c V 2ta )VfJ(Q a + BD c V 2<a ) T ]rl (2.43) 
P =r T (A - Q a V 2 -JC) T P(A - Q a V 2 lC)r 

+ r T (P a + R 2 , a D c C) T R^ a (P a + R 2 a D c C)t (2.44) 

Q =t(A - BR 2 \P a )Q{A - BR^ a P a ) T T r 

+ r(Q a + BDcV 2 ,a)V 2 ~J(Q a + BD c V 2ia ) T T T (2.45) 
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rank P = rank Q = rank QP = n c 


r = (QP)(QP)* 


(2.46) 

(2.47) 


Theorem 2.2 [3]. Suppose D ^ 0 and (A c ,B c ,C c ) is an admissible extremal of the opti- 
mal reduced-order dynamic compensation problem for discrete-time systems. Then, there exist 
nonnegative-definite matrices P,Q,P and Q such that A C ,B C and C c are given by (2.38)-(2.40) 
with D c = 0 for some projective factorization ( G>M,r ) of QP such that conditions (2.42)-(2.47) 
are satisfied with D c — 0. 

Remark 2.3. Theorem 2.1 is a modification of earlier results [2,4]. The primary difference 
is that the P and Q in Theorem 2.1 satisfy the rank conditions (2.46), which parallels the corre- 
sponding continuous-time theory [4,5], whereas the P and Q in [2] and [4] do not satisfy these rank 
conditions. 


The following corollary characterizes the optimal, full-order, discrete-time controller. 


Corollary 2.1. If n c = n x , then one can choose r = r = G = / nx , such that t± = 0 and 
(2.38)-(2.46) reduce to 

A c = A - BC c - B c C - BD c C (2.48) 

Be = Q a V 2 -} - BD C (2.49) 

C c = R 2 ^P* - D c C (2.50) 

D c = R^ a (B T PAQC T + Rj 2 QC T + B T PV l3 )V 2 ~f (2.51) 


where 


P = A t PA + Ri - PjR^ a Pa 

(2.52) 

Q = AQA r + V x - Q.VgQl 

(2.53) 

P — (A — Q a V 2 ~}C) T P(A - QaVfJC) 

+ (Pa + Rl,*D e C) T R 2 \{P a + Rl,aD c C) 

(2.54) 

Q = (A- BR-\P a )Q{A - BR 2 \P a f 

+ (Q a + BD c V 2>a )V 2 -J(Q a + BD c V 2 , a ) T 

(2.55) 

rank P = rank Q = rank QP - n x . 

(2.56) 
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Remark 2.4. Condition (2.56) requires that the LQG controller (A cy B Cy C c ,D c ) have minimal 
order n x . Also, P and Q are not needed to compute the controller but are the closed-loop grammians 
to be used in balanced controller reduction. 

Remark 2.5. Notice that in the full-order case (i.e., n c = n x ), without loss of generality 
one can choose r = G = T = / nx and (2.42) and (2.43) reduce to the standard regulator and 
observer Riccati equations and (2.38)-(2.41) yield the usual LQG expressions. It can be shown 
that (2.44)-(2.46) are equivalent to the requirement that the controller (A c ,R c ,C c ) be minimal. 

Theorem 2.4 [6], Suppose there exists nonnegative definite matrices Q,P, Q and P satisfying 
(2.40)— (2.45) and A Cy B c ,C c and D c satisfy (2.36)-{2.39). Then, the compensator ( A Cy B c ,C c ) is 
an extremal of the optimal fixed-order dynamic compensation problem. Furthermore the following 
are equivalent: 

(i) A is stable 

(ii) (A, V^) is stabilizable 

(iii) (A, iH) is detectable. 

In addition, 

(A c , B c ) is controllable <=> A c + B c CG T is stable (2.57) 

(A c ,C c ) is observable <*=> A c + rBC c is stable. (2.58) 

In the homotopy algorithms to be subsequently defined the optimal projection equations (2.42)- 
(2.45) due to their relationship to standard LQG equations can be used to give insights into the 
development of initializing controllers. However, the homotopy algorithms will be based directly 
on the gradient of the cost functional. 
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3. Review of Homotopy Methods 

A “homotopy” is a continuous deformation of one function into another. Over the past several 
years, homotopy or continuation methods (whose mathematical basis is algebraic topology and 
differential topology [1]) have received significant attention in the mathematics literature and have 
been applied successfully to several important problems [2-7]. Recently, the engineering literature 
has also begun to recognize the utility of these methods for engineering applications (see e.g. [8- 
10]). The purpose of this section is to provide a very brief description of homotopy methods for 


finding the solutions of nonlinear algebraic equations. 

The reader is referred to [7,8,11,12] for additional details. 

The basic problem is as follows. Given set 0 and $ contained in IR n and a mapping F : 0 — ► 
find solutions to 

F{6) = 0. (3.1) 

Homotopy methods embed the problem (3.1) in a larger problem. In particular let H : 0 x [0, 1] — ► 
IR n be such that: 

1 ) H(e,l) = F{9). (3.2) 

2) There exists at least one known do € IR n which is a solution to H(-,0) = 0, i.e., 

H{9 o,0) = 0. (3.3) 

3) There exists a continuous curve (0(A), A) in IR n X [0, 1] such that 

H{6{ A), A) = 0 for A € [0,1] (3.4) 

with 

(*(0),0) = (fl o ,0). (3.5) 


4) The space 0 x [0, 1] has a differential structure so that the curve (0(A), A) is differentiable. 

A homotopy algorithm then constructs a procedure to compute the actual curve a such that the 
initial solution 0(0) is transformed to a desired solution 0(1) satisfying 

O=tf(0(l),l) = F(6l(l)). (3.6) 


Differentiating if (0(A), A) = 0 with respect to A to obtain Davidenko’s differential equation 



OH dd dH 
~d9dX + d\ ~ ' 

(3.7) 
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Together with 0(0) = 0 O , (3.7) defines an initial value problem which by numerical integration from 

0 to 1 yields the desired solution u(l). Some numerical integration schemes are described in [11,12]. 
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4- The Homotopy Map and It’s Jacobian 
If we define 


0± 


(4.1) 


[■ vec(A c ) ' 
ve c(B c ) 
vec(C c ) 

Lvec(D c )J 

then the cost functional of Section 2 can be expressed as J(0). The homotopy defined in this 


section is based on finding 6 satisfying 


o = m = 


It is useful to recognize that 


vJW T * % = 


yec dA c 


vec 

vec 


dj_ 

dB c 

dj_ 

dC c 


vec ^z_ J 

yec dD c J 


(4.2) 


(4.3) 


Expressions for the partial derivatives > "§b~ c > » an d 1)D^ are derived in Appendix A. Here, 
we cite only the final results. First, we assume that P, Q and Z satisfy 


P = A 1 PA + R 
Q = AQA t + V 
z = qa t p 

and note that P,Q and Z have the partitioned forms 


P = 


m 

p t 

M2 


P 12 
P 22 


Q = 


Qu Q 12 
Qn Q 22 


z = 


Z_w 
Z 21 


Z 12 

Z22 


(4.4) 

(4.5) 

(4.6) 

(4.7) 


where the (1.1) and (2.2) blocks of each matrix are respectively n x n and n c x n c . With this in 
mind, the cost derivatives are given by 


M. = 2Z? 
»A U 


22 


dj 


= 2(P£V„ - P&BDcVi + P2iB c V 2 
+ zJ 2 c T - zJ 2 cJd t ) 

= 2 ( — R.J 2 Ql2 + RiD'CQu + RlC c Ql2 
4- B t Z% - D r BjzJ 2 ) 

= 2 (-Rj 2 QuC T + R 2 D c CQ u C t + R 2 C c Qj 2 C T 

(j D Q 

- B J P u Vn + B T PnBD c V 2 - B r PuB c V 2 


dB c 

dj 

dC c 


(4.8) 

(4.9) 
(4.10) 


- B t Z^C T + R 2 D c V 2 ). 


(4.11) 
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Definition of the homotopy map h(ff, A) 


To define the homotopy map we assume that the plant matrices (A, 0, C, D), the cost weighting 
matrix (R\, R 2 ,Ru) and the disturbance matrices (Vi, V 2 , Vi 2 ) are functions of the Homotopy 
parameter A £ [0, 1]. In particular, it is assumed that 


A(A) S(A)1 
C( A) D(A)J 


* £]+*([* B ’ 


Co 


Cf Df 


Aq Bq 
C o Bq 


)• 


(4.12) 


where 


and Lx# and Lrj satisfy 


where 


and Lv t o and Lvj satisfy 


MX)} “ L * WL *M 

Lr( A) = LR t o + A (i/R,/ — LRfl) 


Lr$Lr 0 — 

LrjLrj = 


Ri,o Ru,o 

1^12,0 ^2,0 J 

R\,J R\2J 
f?2 ,/ 


L-^12./ 


F,(A) V 12 (A) 

^12 (A) V 2 T (A) 


= L v (X)Ll(X) 


Lv{ A) = Lv,o + A (Lvj - Lv,o ) 


I'Ve Lj / 0 


[ Vlfl Vl2,0 

V$, a ^.oj 

T rT _ f V,., V UJ 

Lv ’ L VJ -\v&j V,.o 


(4.13a) 

(4.13f») 

(4.13c) 

(4.13d) 

(4.14a) 

(4.146) 

(4.14c) 

(4.14</) 


Note that (4.12)-(4.14) imply that >1(0) = Ao and j4(1) = A/, B( 0) = Bo and .0(1) = 0/, etc ... 
and it is understood that Aj,Bj,... were referred to previously simply as A, 0, — The change in 
notation is simply for convenience. 


The homotopy map h(6 , A) is defined by 


h(0, A) 


'vec (H Ac ( 6 , A))' 

vec(0 Bc (0,A)) 

vec(0 Cc (0,A)) 

_vec(0 Dc (0,A)). 


(4.15) 
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where 


^ c ( 0,A) = 2Z 2 t 2 (4.16) 

H Bc ( 9 ,X) = 2 (P&V12 - P? 2 bd c v 2 + P 22 B c v 2 

+ Zj 2 C J - Zj 2 CjD T (4.17) 

tfc.(M) = 2(— + RiD c CQu + RiCcQh 

- B^ Zj 2 - D T Bj zj 2 (4.18) 

= 2 (— 

- B t P„V 12 + P t P„£I> c V 2 - B T Pi2PcV 2 

-P t Z 1 t 1 C t + P2PcV 2 ). (4.19) 


dh dh 
d9 dA 


The Jacobian of the homotopy map 

We now consider that computation of ^h(6, A)^, the Jacobian of h{9, A). Note that 

v/i(M) T = 

Recalling that 6 is defined by (4.1), such that for some integers k and t , 0 3 is given by 

Oj = &c,kt'> @3 = — C c,kfi or — d c% kf 

It follows from (4.6) that is of the form 

• * • ve <aZT t HB ‘ ) • ' • vec ( ) • ' • vec ( : Hb ^- ■ ■ ve <a£j ; HB ‘ ) 

• • • vec (a^7^ ) • • • vec (WT t Hc ‘ ) • • ' vec ( a^I7^ c ‘ ) • ■ • vec (^7^ ) 

...vec(^P 0 J..- v ec(^^ e )...vec(^ 7 P D J...vec( 3 ^ 7 P Dc )J 


dh 

d0 


and can be expressed as 


dh 

d\ 


vec(£H Ac ) 
vec(^P B J 
vec (jiH Cc ) 
vec (jx H d c ) 


(4.20) 


(4-21) 


(4.22) 


Below, we develop explicit expression for the derivative terms appearing on the right hand 
(4.22) and (4.23). We use the notation 

M U) £ ™ 

d9j 

m* 9M 

M= dA* 


(4.23) 
sides of 

(4.24) 
(4-25) 
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Differentiating (4.4)-(4.6) with respect to 0j yields 


P U) = A T P U) A + {A U)T PA -1- A t PA U) + fl (j) ) 

(4.26) 

QU) = AQ^A t + (A<»QA T + AQA^ T + V^) 

(4.27) 

Z (j) = Q<»A t P -1- QA (i)T P + QA T pM 

(4.28) 

where expressions for the derivatives A^\R^ and are given by (B.20)-(B.28) of Appendix B. 

Similarly, differentiating (4.4)-(4.6) with respect to A yields 

•_T * * 

P = A 1 PA + (A PA + A T PA + R) 

(4.29) 

Q - AQA T + (AQA t + AQA + V) 

(4.30) 

Z = $a t p + qa P + qaP 

(4.31) 


where expressions for A y R and V are given by (B.29)-(B.33) of Appendix B. 

Before presenting the desired derivative expressions we define 

H' Ac (Z^) = 2Zif (4.32) 

H' Bc (pl j Kz^) = 2{P[fv n - P[f BD C V 2 + P^B C V 2 

+ Z[fc T - Z { 2 i )T CjD T ) (4.33) 

H' Cc {Q (i \Z^) = 2(-Rj 2 Q$ + R 2 D c CQ\ '$ + R 2 C c Q$ 

- B T Z ( 2 j )T - D T BjZ<i )T ) (4.34) 

H' Dc (P (i) ,Q U) ,Z<») = 2(-Rj 2 QifC T + R 2 D c CQ^C t + R 2 C c Q\fc r ) 

- B T Z\i )T C T ) (4.35) 


Notice that the right hand sides of (4.32)-(4.35) are essentially identical in form to the right 
hand sides of (4.16)-(4.19). The difference is that P,Q , and Z have been replaced by P^\qU) 
and Z ^ and the last term ( 2R 2 D C V 2 ) in (4.19) has no counterpart in (4.35). 


ivatives with Respect to b. 


Differentiating (4.16)-(4.19) with respect to b c< kt(= Oj) gives the following. 


dHA< 

db Ci kt 

dHjh 

db cM 


H' Ac {Z U) ) 

H' Bc (P^\Z^) + 2P 22 E[ k /J n V 2 


(4.36) 

(4.37) 
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|^- = H' Ct (Ql»,Zl») - 2D T Ei t * ) nc Z£ 

O&cM 

= H' Dc (P U \Q U \Z U) ) ~ 2B 7 P u E^ n V 2 . 

oB Cy kt 


Derivatives with Respect to w 

Differentiating (4.7)-(4.9) with respect to 0j(= c cM ) gives the following. 

= H' Ac (Z^) 

VC Ct ki 

= H' b (P u \Z (j) ) - 2 Zj 2 E { n ‘' k) n D T - 2(fi (i) Zfi + SlZSl ij) )B c 

dc c ,kt 


dH Cc 

dc cM 


H'cSQ u \Z^) + 2R 2 El k * nc Q 22 


dHp £ 
dc c ,kt 


H' Dt {P (i \Q U) ,Z U) ) + 2R 2 E ( n k i ) nc Qj 2 C 


Derivatives with respect to d c 1 j_ 


Differentiating (4.7)-(4.9) with respect to d c<k i gives the following. 


dH± 

dd Cj ke 

dH^ 

dd Cf kt 

dHc, 

dd^kt 

OHl > c 

&d Ct k£ 


H’ Ac (Z »>) 

H' Bi {P^\Z^) - 2P? 2 BE (k y n V 2 
H' Cc (Q^,Z {j) ) + 2R 2 E { n k ^ n CQ 12 

H' Dc (P {j \Q U) ,ZM) + 2{R 2 E (k y n CQ u C T + B 7 P n BE (k ^ n V 2 

+ R.E^ilv,). 


(4.38) 

(4.39) 


(4.40) 

(4.41) 

(4.42) 

(4.43) 


(4.44) 

(4.45) 

(4.46) 

(4.47) 


Derivatives with Respect to A 
Differentiating (4.7)-(4.9) with respect to A gives 


dH Ac 

ax 

dH Bc 

dX 



(Z) 


= h' Bc {P,'z ) 

+ 2(P 2 t 1 Vi 2 - P 7 2 bd c v 2 - P? 2 BD C V 2 + P 22 B c V 2 ) 


dH Cc 

dX 


+ z 7 2 c 7 - z 7 2 cji) 7 ) 

H'cSOj) 


(4.48) 


(4.49) 
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+ 2 ( — RJ2Q12 + R-2DcCQl2 + RiDcCQu + R-2CcQi2 
+ i? T z£ - D t BJzJ 2 ) 


dHp c 

3A 


= H' Dc (P,Q,Z) 


+ 2(-j?7 2 g„c T - j&QhC t 

+ R 2 D c CQi\C T + ^■^cC'QnC' 1 ' + i?2 D c CQ n 

+ + r 2 c c qJ 2 c t ) 

+ 2(-B T P n V u - £ T P n Vi2 


+ B T P n BD c V 2 + B T PnBD c V 2 + B r P u BD c V 2 
- B t P 12 B c V 2 - B T P n B c V 2 ) 

- 2B r Zj 1 C T - B T Z n C T 
+ 2(R 2 D C V 2 + * 2 ^ 2 ) 


(4.50) 


(4.51) 


where from (4.12)-(4.14) 


where 


where 


'A B 


Af — Ao Bf — Bo 

C D 


Cj — Co Cj — Co 



•fti2 

#2 


= LrLr 


Lr - LRJ - IyRfl 



Vl 2 

V 2 


= LyLv 


Lv = Lvj - Lv,o- 


(4.52) 

(4.53a) 

(4.536) 

(4.54a) 

(4.546) 


The homotopy Jacobian can now be computed using (4.20) with (4.22) and (4.23) and (4.32)- 
(4.51). Note that the primary computations involve the computation of the solutions of the Lya- 
punov equations (4.26), (4.27), (4.29) and (4.30). Significant computational savings can be made 
by solving these Lyapunov equations in a basis in which the A matrix is diagonal (or nearly diago- 
nal). This requires transforming the corresponding forcing terms into this basis. But it is seen by 
(B.20)-(B.33) of Appendix B that these forcing terms are very sparse. Hence this transformation 
does not have to be expensive. In addition, it is required that the solutions of the Lyapunov equa- 
tions be transformed into their original basis before substituting into the expressions (4.32)-(4.51). 
A close examination of these expressions shows that for problems in which n y << n x ,n u << n x 
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and/or n c « n x significant computational savings can be made by not actually performing the 
matrix multiplies to transform the solutions into their original basis until after substituting the 
transformations into (4.32)-(4.51). Appendix H gives the details of efficient computation of H$ for 
the corresponding continuous- time problem. A nearly identical procedure has been implemented 
for the discrete-time problem considered in this report. 
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5. Reduction of the Dimension of the Controller Parameter Vector ( 9 ) 

The homotopy function H(9, A), described earlier, was defined to solve the H 2 optimal reduced- 
order dynamic compensation problem for discrete-time systems. The vector 9 was defined such 
that it contained each of the elements of the controller matrices, A c , B c and C c . However, for 
computational efficiency it is desired that 9 be as small as possible. Hence, we desire to represent 
the controller matrix with the fewest parameters possible (i.e., we desire 9 to have the smallest 
dimension possible). The minimal number of parameters p m i n with which a compensator can be 
represented is given by [1,2] 

Prnin = n c ( n « + n v) ( 51 ) 


One canonical form which allows representation of a controller with a minimal number of 
parameters is the modal form described in [3]. This form will be called here the Second-Order 
Polynomial (SP) form. For this parameterization a triple (A c , B c ,C c ) has the following structure. 

A c = block- diag{j4 c ,i,;4 c ,2 . . .,4r) (5-2) 


where A c ,i G IR 2 * 2 for i G {l,2,...,r} and each A Cii (with the exception of A C<T if the row 
dimension of A c is odd) has the form 


Ac t i — 


0 

c, X 


1 


(5.3) 


to allow for either a complex conjugate set of poles or two real poles. B c is completely full and 


C c — [Cc,l, C C) 2, • • • ,C Cj r] 


(5.4) 


where C Ci , has the form 


C 


c,r — 


‘1 

* 


* 




* 


(5.5) 


The controller canonical form described in [4,5] also allows representation of a controller with 
a minimal number of parameters. For single-input, single-output (SISO) systems in controller 
canonical form the A c matrix is a companion matrix. In particular, A c has the form 


A c = 


-0100 ••• 0 ' 

0 0 10 0 

0 0 0 1 0 . 


* * * * • • • * J 


(5.6) 
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In addition, 



0 

1 


(5.7) 


and C c is completely full. A dual form of the controller canonical form is the observable canonical 
form [5]. 


It is also possible to represent the controller in a basis where the number of free parameters p 
satisfies 

Pmin < p < Pmax = «c(«c + + n y ). (5.8) 


One such basis is the tridi agonal basis [7-11] in which the controller state matrix is constrained to 
have nonzero elements only on the diagonal, the super-diagonal and the sub-diagonal. That is, 


A c = 


* 

* * 

* * 



(5.9) 


L * *J 

B c and C c are completely full. For this form the number of free parameters is given by 


P — Pmin 4" (3n c 2) 


A common feature of each of the above bases is that they are described by simply constraining 
certain elements of the controller (or plant) matrices to constant values (e.g., 1 or 0) while allowing 
the remaining parameters to have arbitrary values (A c , B cy C c ). Hence, the corresponding parameter 
vector (0p) y gradient vector (Je.p) and Hessian matrix (H$ tP ) are given by 


q* 

Ph 

II 

(5.10) 

Je,p — r Je 

(5.11) 

H», p = TH e T r 

(5.12) 


where T is an elemental matrix (i.e., each row has only one nonzero element and this element has 
unity value). It should be noted here that He tP can be computed more efficiently than shown in 
(4.64). Since it is not necessary to construct the large Hessian H$ to compute the smaller Hessian 
He lP . 
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6, Overview of the Homotopy Algorithm 

This section describes the general logic and features of the homotopy algorithm for H 2 opti- 
mal reduced-order control. It is assumed that the designer has supplied a set of system matrices, 
5/ = describing the optimization problem whose so- 

lution is desired. In addition, it is assumed that the designer has chosen an initial set of related 
system matrices So = {Ao,Bo,Co,Do,Rifi, -# 2 , 0 , Fi.o, I^.o, Fn.o) that has an easily obtained opti- 
mal controller (A c ,o, -Sc.o, C C fl,D c ,o) of order n c . The initial system So can be chosen to correspond 
to a low-authority control problem as described in Appendix I since if Ru, or V% 0 are of the appro- 
priate structure and Aq is stable, the corresponding LQG controller is nearly nonminimal and can 
hence be reduced to a nearly optimal n‘ h order compensator using, for example, balanced controller 
reduction [1]. The reader is referred to Appendix I for additional details. 

Below, we present an outline of the homotopy algorithm. This algorithm describes a predic- 
tor/corrector numerical integration scheme. There are several options to be chosen initially. These 
options are enumerated before presenting the actual algorithm. Notice that each option corresponds 
to a particular flag being assigned some integer value. 

Controller Basis Options : 

basis = 0. No basis (i.e., all elements of the controller matrices are considered free.) 

basis = 1. Tridiagonal Basis. 

basis = 2. Second-Order Polynomial Form. 

basis = 3. Controller Canonical Form. 

Note that for basis = 0 or 1, p > p min while for basis = 2 or 3, p = p m i„. 

Prediction Scheme Options : 

Here we use the notation that Ao,A_i, and Aj represent the values of A at respectively the 
current point on the homotopy curve, the previous point and the next point. Also, 6 p ' = dO p /d\ 
and is the solution of Davidenko’s differential equation (4.7), rewritten here as 

H e , p 0’ p ( A) + H X = 0. (6.1) 

If p = p m in, He, p is generally invertible, then 0p(A) is given exactly by 

0' p (X ) = (6.2) 
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If p > p m i n , then Hb, p generally has rank p m ;„ and 0J,(A) is approximated by the least squares 
solution of (6.2) or 

^ =v [T l] uT (6 - 3) 

where it is assumed the Hg tP has the singular value decomposition 

H,, p = U JJf t , (6.4) 

Note that for p = p m i n (6.3) and (6.4) are equivalent. 

pred ss 0. No prediction. This option assumes that 0 p (Ai) = 9 p ( Ao). 

pred = 1. Linear prediction. This option assumes predicts 0p(Aj) using only 0 P (A O ) and 0 p '(Ao). 
In particular, 

9 P ( Ai) = ^p(Ao) + (Ai — Aq)0 p '(Ao) (6-5) 

pred = 2. Cubic spline prediction. This option predicts 0 p (Ai) using 0 p (A o ), 9 p (\o), 0p(A_i) 
and ^(A^i). In particular, 

0p(Aj) = gq + aiAi + U2A1 2 + U3A1 3 (6-6) 

where ao, a\ y a 2 and a 3 are computed by solving 

'1 0 1 0 1 r0p(A-i) 

r i A_i 1 Ao 1 _ 5'(A_i) 

[a 0 01 a 2 asj ^ 2Aj ^ ^ 

lAi, 3A1 x Ag 3AgJ L^(Ao) 

Note that if this option is chosen, then at the initial algorithm prediction step 0 p (A_i) and 
0p(A_i) are not available, in which case linear prediction is used. 



Correction Options: 


Here we assume that the homotopy parameter has a fixed value Ao- The vector 9 V represents 
the current approximation of the parameter vector at A = A<>. Each of the options corresponds to 
updating 9 P using the formula 

9 P ^6 P + A 9 P (6.8) 


where 


A0 P — 


(6.9) 


for some choice of G$ tP . 
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cott = 1. Newton correction. In this option, if p = Pmim 

G e>P = H e , p ~ x (6.10) 

while if p > Pmin i 

G e , p = V(£ 2 +a 2 /) -1 Ef/ T (6.11) 

where a is some (small) scalar and (I/, V, E) denote the singular value decomposition of H t, p 
such that 

H e<p = UZV T . (6.12) 

It can be shown that if Ge, P is given by (6.11), then A 0 P minimizes 

+ 6 p || 2 + a 2 ||A5 p || 2 ]. (6.13) 

Hence, A 0 p is essentially a “Newton correction” that is relatively insensitive to singularities or 
near singularities in the Hessian, H$ <p . 

corr = 2. Quasi-Newton correction. In this option, Ge, P denotes the estimate of * using 
only gradient and cost information. For the algorithm presented here the BFGS inverse Hessian 
update is used [2]. 

Outline of the Homotopy Algorithm 

Step 1. If basis > 1, then transform the initial controller (A Ci o,Bc,o>C'c,o) to the chosen basis 
and let 0 o ,p be the corresponding vector of free parameters. 

Step 2. Initialize loop = 0, A = 0, AA € (0,1], S = 5 0 , B p = 0 O , P and compute the cost J 
and the cost gradient J e, P corresponding to S and the controller described by 6 P . 

Step 3. Let loop = loop+1. If loop = 1, then go to Step 5. Else, continue. 

Step 4. Advance the homotopy parameter and predict the corresponding parameter vector 6 
as follows. 

4a. Let Ao = A 

4b. Let A = Aq 4- AA. 

4c. If pred > 1, then compute 0p(A o ). 

4d. Predict # P (A) by using the option defined by pred. 
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4e. If the normalized gradient J 9 tP \\G » tP \\j ^0 P satisfies some preassigned tolerance, then 
continue. Else, reduce AA and go to Step 4b. 

Step 5. Correct the current approximation 8 p to the optimization problem defined by S using 
the option defined by corr until the normalized gradient, 

(6.14) 

A Op 

satisfies some preassigned tolerance. 

Step 6. If A = 1, then stop. Else, go to Step 3. 

The above algorithm assumes monotonicity of the solution curve as a function of the homotopy 
parameter A. However, it is not difficult to modify the algorithm so that the variable parameter is 
the arc length as discussed in [3,4] since this modification would still only require the computation 
of H$ and H\. The modified algorithm would not require monotonicity of the solution curve. 
However, so far in our computational experience the solution curve has always been monotonic. 

Note that if p = p m ; n and corr = 1, then the corrections of Step 5 correspond to Newton 
corrections. Hence if the prediction tolerance used in Step 4 is sufficiently small, then, entering 
Step 5, 6 p will be close enough to the optimal value 0* so that the quadratic convergence properties 
of Newton’s method [2] can be realized. In practice, this quadratic convergence property is not 
always realized due to numerical ill- conditioning associated with the minimal parameterization of 
the controller. This ill-conditioning is illustrated and discussed further (in the context of continuous- 
time systems) in Appendix H. 
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7* A Design Example and Some Rules of Thumb 

This section illustrates the design of a reduced order compensator for an axial beam with 
four disks attached as shown in Figure 7.1. This example has been considered in several previous 
publications [1-7] and has become a standard benchmark example. The section closes with some 
general rules of thumb that will aid the control designer in most efficiently utilizing this algorithm. 

The basic control objective for the four disk problem is to control the angular displacement 
at the location of disk 1 using a torque input at the location of disk 3. It is also assumed that a 
torque disturbance enters the system at the location of disk 3. An 8th order discretized model of 
the fourdisk plant with nominal performance weights and disturbance covariances is generated by 
the function diskmod. 

The design philosophy adopted here is that the scaling q2 of the nominal control weight R 2 
and the nominal sensor noise covariance V 2 are simply design knobs used to determine the con- 
troller authority. The system costs are computed assuming that V2 = 0. This general philosophy is 
actually motivated by insights into LQG theory. However, it will suffice here to simply note that 
this philosophy was used successfully on two hardware experiments involving control design and 
implementation [10-13]. 

It is now assumed that we are in the MATLAB environment. In what follows the reader is 
walked through the design process for a 4th order controller. The command sequences are presented 
after the prompt and after the commands some of the resultant output is displayed. Explana- 
tory text is interspersed to clarify the motivation of the command sequences and the interpretation 
of some of the output. 

We begin by using diskmod to generate the design plant and nominal weights and covariance, 
diskmod 

discretizing a, b, and vl 
The following variables are now loaded into memory. 
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>> who 

Your variables are: 

a c rl r2 vl v2 

b d rl2 ts vl2 

We choose to display numerical data using the following format. 

» format short e 

We now begin the search for an authority level that will give us a nearly optimal controller by 
balanced controller reduction. We commence this process by choosing the initial scaling q20 of R 2 
and V 2 as follows. 

» q20 « .1; 

We use dlqg to design an LQG controller and then check the eigenvalues of the product phat*qhat 
to see if their is any gap between the 4th and 5th eigenvalues (ordered in descending order of 
magnitude). Note that the warning after the call to dlqg in this case is not important. 
[ac,bc,cc,dc,costslqg f phat ,qhat] = ... 

dlqd(a,b,c ,d,rl ,q20*r2,rl2,vl ,q20*v2,v!2, 1) ; 

Warning: Q is not symmetric and positive semi-definite 

> -sort (-eig(phat*qhat) ) 
ans * 

1.3554e+01 
1.2377e+00 
8 . 0067e-01 
1 . 3682e-01 
1.0451e-01 
1.7751e-02 
1.0585e-02 
4.9872e-03 

Note that their is no gap between the 4th and 5th eigenvalues indicating that balanced controller 
reduction will probably not yield a nearly optimal reduced-order controller. However to verify this 
we will actually compute a 4th order controller using balanced controller reduction and compare 
it’s cost with that of the LQG compensator, which is contained in the vector costslqg. 

> [ac,bc,cc,dc] = baler ed(ac, be, cc,dc,phat ,qhat ,4) ; 
>costs=dqcosts(ac,bc,cc,dc,a,b,c,d,rl ,q20*r2,r!2. 
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vl,q20*v2,vl2); 

> costslqg 
costslqg = 

2.1047e-02 1.6027e-02 0 4.2135e-04 3.7495e-02 

5 > costs 
costs * 

2.5838e-02 2.06230-02 0 4.2135e-04 4.6882e-02 

The total cost of the LQG compensator is 3.7495e-02 while the cost of the reduced order controller 
is 4.6882e-02. The vast differences in these costs is another indication that the reduced order 
controller is not nearly optimal. We will now repeat the above process with a higher value of q20, 
i.e. the LQG controller is of lower authority. 

> q20 = 10; 

>>[ac,bc,cc,dc # costslqg,phat ,qhat] = ... 

dlqg(a ( b ( c,d,rl ,q20*r2,rl2 , vl ,q20*v2 ( vl2, 1) ; 

Warning: Q is not symmetric and positive semi-definite 

> -sort (-eig(phat*qhat) ) 
ans = 

4.1835e+02 
2 . 1594e+00 
5 . 6033e-01 
4 . 3915e-01 
4.6232e-02 
3.7616e-02 
4.4658e-03 
4.4134e-03 

> [ac,bc,cc,dc] = balcred(ac,bc,cc,dc,phat ,qhat ,4) ; 
^costs^dqcostsCac^c^c^Gja.b^c.d.rl ,q20*r2 ,rl2 , vl , 

q20*v2 , vl2) ; 

> costslqg 
costslqg = 

2 . 0046e-01 4.1697e-01 0 2.7474e-03 6.2018e-01 

>> costs 
costs = 
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2.0137e-01 4.1748e-01 0 2.7474O-03 6.21600-01 

This time there is an order of magnitude gap between the 4th and 5th eigenvalues of phat*qhat 
and the costs of the LQG and reduced-order controllers are nearly identical. This indicates that 
the 4th order balanced controller is nearly optimal. This deduction could also be made by gener- 
ating a performance curve for the LQG controller (by varying q20) and superimposing it with the 
performance curve for the corresponding 4th order balanced controllers as shown in Figure 7.2. If 
for a given q20 the two controllers have essentially the same state and actuation costs then the 4th 
order balanced controller is probably nearly optimal. 

If the 4th order balanced controller is nearly optimal then by using a few Newton corrections 
(say, 1 to 6) we should be able to converge to the optimal controller (practically, the controller 
that satisfies a small tolerance on the normalized norm of the cost gradient). This is verified 
below. Function nwtpar is used to intialize the algorithm parameters to their default values while 
nwtprint is used to display these default parameters. 

>> par = nwtpar; 

> nwtprint (par) ; 

1. Will print intermediate results. 

2. gradient prediction tolerance = 1.00000e-05 

3. gradient correction tolerance = 2.00000e-08 

4. gradient final tolerance = 2.00000e-08 

5. minimum homotopy step size * 1.00000e-06 

6. maximum number of corrections allowed « 10.000000 

7. Will use Hessian from last correction for prediction. 

8. Will not use line search. 

9. Will let program run automatically. 

10. initial step size * 1.000000 

11. No basis is assumed for the controller. 

At this time the user has the option of changing any of the default parameters. However, we will 
be content with them. The default parameters will also be printed by dnwthom. In the following 
call to dnwthom the initial and final system parameters are identical so that the algorithm will 
only perform correction loops. 

» [ac ,bc,cc,dc,val ,par] = dnwthom(ac,bc,cc,dc, ... 
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a,b,c,d,rl f q20*r2,rl2,vl ,q20*v2,v!2, . . . 
a,b ,c,d,rl ,q20*r2 ,rl2 , vl f q20*v2, vl2 ,par) ; 

1* Will print intermediate results. 

2. gradient prediction tolerance « 1.00000e-05 

3. gradient correction tolerance * 2.00000e-08 

4. gradient final tolerance = 2.00000e-08 

5. minimum homotopy step size s 1.00000e-06 

6. maximum number of corrections alloved = 10.000000 

7. Will use Hessian from last correction for prediction. 

8. Will not use line search. 

9. Will let program run automatically. 

10. initial step size = 1.000000 

11. No basis is assumed for the controller. 

Computing Initial Hessian. . . 

Inverting Hessian... 

***** INITIAL PARAMETERS ***** 

normalized normalized normalized 
cost costO-cost gradient gradient del-theta 

6 . 21598e-01 0 . 00000e+00 7 . 35784e-07 5 . 54898e-010 . 00000e+00 

The algorithm is still in process but we note here that the initial 
normalized gradient value of 7.36e-07 is fairly small. The general 
rule is that values < to about 2.0e-08 are very 
small . 

********** CORRECTING ********** 

**** lambda = 1.0000e+00 **** 

Correction Iteration 1 

Computing Hessian. . . 

Inverting Hessian. . . 

** correcting: i = 1.000000, lambda = 1.0000e+00 ** 

normalized normalized normalized 

cost costO-cost gradient gradient del-theta 

6 . 21598e-01 0.00000e+00 7.35784e-07 5.54898e-01 0.00000e+00 

6 . 21515e-01 1.33389e-04 7.34790e-08 5.S4259e-02 1.92813e-03 
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The normalized gradient correction tolerance is: 2.00000e-08 

With the algorithm still in progress we note that the top line denotes the initial convergence 
parameters before the first correction while the bottom line denotes the value of the convergence 
parameters after the correction. It is seen that both the cost and gradient were improved (i.e., 
decreased by the first correction iteration). However the normalized gradient is still not below its 
maximum tolerance of 2.0e-8. 

Correction Iteration 2 

Computing Hessian... 

Inverting Hessian... 


** correcting: i = 2.000000, lambda = 1.0000e+00 ** 

normalized normalized normalized 

cost costO-cost gradient gradientdel-theta 


6.21515e-01 1.33389e-04 7.34790e-08 5.54259e-02 1.92813e-03 

6.21515e-01 5.92847e-07 2.33195e-10 1.77478e~04 1.86134e-04 

The normalized gradient correction tolerance is: 2.00000e-08 

doubling step size to 0 

**** Exiting DNWTHOM with lambda=l. **** 

The correction of the initial 4th order balanced controller converged in 2 iterations. This clearly 
indicates that the balanced controller was nearly optimal. The controller (ac,bc,cc,dc) is now the 
optimal 4th order controller for the scale factor q20. 

We now set q2=l and use the dnwthom to deform the initial controller into a higher authority 
controller. We show some of the beginning output of dnwthom. 

[ac,bc,cc,dc,val,par] * dnwthom(ac,bc,cc,dc, ... 

a,b,c,d,rl ,q20*r2,rl2,vl,q20*v2,v!2, . . . 
a,b,c, d,rl ,q20*r2,rl2,vl ,q20*v2, v!2,par) ; 


***** INITIAL PARAMETERS ***** 

normalized normalized 
cost costO-cost gradient gradient 


normalized 

del-theta 


6.21515e-01 O.OOOOOe+OO 2.58673e-10 1.73424e-04 0.00000e+00 
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********** PREDICTING ********** 

**** lambda = 0.0000e+00 **** 

** dlambda * 1.0000e+00 ** 

Computing Pseudo-Inverse of Hessian! 
number of sing. vals. retained s 9.000000 

** predicting: lambda * 1.00000e+00, dlambda = 1.00000e+00, ** 

normalized normalized normalized 

cost costO-cost gradient gradient del-theta 

6 . 21515e-01 0.00000e+00 2.58673e-10 1.73424e-04 0.00000e+00 

1.72524e-01 2.60247e+00 8.79235e-05 5.89473e+01 1.95698e-02 

The normalized gradient prediction tolerance is: 1.00000e-05 

!! adjusting step size!! 

dlambda = 5.0000e-01 


** predicting: lambda = 1.25000e-01, dlambda = 1.25000e-01, ** 

normalized normalized normalized 

cost costO-cost gradient gradient del-theta 

4 . 79497e-01 2.96182e-01 1.46945e-05 9.85177e+00 4.89245e-03 

5.48829e-01 1.32438e-01 4.47514e-06 3.00030e*00 2.44622e-03 

The normalized gradient prediction tolerance is: 1.00000e-05 

********** CORRECTING ********** 

**** lambda = 1.2500e-01 **** 

Correction Iteration 1 

Computing Hessian... 

Inverting Hessian... 

condH(l) * 3.36501e+07 [Hessian condition number] 

condH(2) = 9.61511e-01 [free parameter singularity] 

condH(3) = 2.05027e+00 [dthetap ratio] 

** correcting: i = 1.000000, lambda = 1.2500e-01 ** 

normalized normalized normalized 

cost costO-cost gradient gradient del-theta 


5.48829e-01 0.00000e+00 4.47514e-06 3.00030e+00 O.OOOOOe+OG 

5.48807e-01 4.01362e-05 1 . 10053e-07 5.74753e-02 1.54763e-04 
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**** Exiting DNWTHOM with lambda=l. **** 

We now use valprint to check the performance parameters recorded in the vector val. 
>valprint(val) 

1. final homotopy parameter value * 1.00000e+00 

2. total # of Hessian calculations * 33.000000 

3. minimum # of corrections for fixed lambda = 1.000000 

4. maximum # of corrections for fixed lambda * 3.000000 

5. minimum homotopy step size - 3.12500e-02 

6. maximum homotopy step size = 1.25000e-01 

The costs plotted in Figure 7.2 are computed using dqcosts as follows. Note that in the input 
arguments v2 is set to zero. 

>costs = dqcosts (ac, be, cc,dc,a,b,c,d,rl ,r2,r!2,vl ,0,vl2) ; 

^>costs 


costs 

= 





5.5053e-02 

7.8950e-02 

0 

0 1.34000-01 

The optimal controller is listed below. 



> ac 





ac = 

9 . 6632e-01 

4.67900-02 

-1.19160-02 

-5.39260-03 


-3 .0758e-02 

9.60530-01 

7.4261O-03 

-4.10060-04 


2 . 9294e-03 

-8.69590-03 

9.93350-01 

8.7504O-02 


1.3349e-03 

5.57980-03 

-8.8247e-02 

9.89800-01 

> be 





be = 

-2.9757e-02 
9.2577e-02 
-3 . 3036e-02 
-2.8086e-02 




> cc 





cc = 

-8.6600e-02 

6.44350-02 

1.69830-02 

-2.97950-02 
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> dc 
dc = 

3.4606e-02 

Some Rules of Thumb 

1. Choose the initial weights (Ri t o*R 2 ,(hRi 2 t o ? V| t o, V^o, Vi2,o) an d the final weights (Rij* ^2,/> 
Ri 2 j,VijyV 2 j>Vi 2 ,f) so that along the homotopy path the regulator and estimator poles have 
the same order magnitude. That is avoid situations where the estimator poles are very fast while 
the regulator poles are slow or vice versa. The algorithm will converge in these cases but the 
convergence tends to be slow. 

2. Our experience indicates that no constraints on the controller basis appear to yield better 
numerical robustness than constraining the basis to tridiagonal form or some other basis. When 
attempting a 6th order controller for the four disk problem constraining the controller basis to 
tridiagonal form yielded very poor numerical robustness. However, when the controller basis was 
left unconstrained the algorithm performed very well. This phenomena is discussed further in 
Appendix H. 

3. For better control don’t take huge steps between the initial and final system parameters. For 
example don’t try to go from very low control authority to very high control authority all at once. 
Take “reasonable size” increments. You may want to adjust the tolerances along the way to increase 
the algorithm efficiency. 
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8. HADOC Toolbox Reference 

balcred 8-6 

beam 8-7 

bodeplot 8-8 

clp,dq,clz 8-11 

dclp,dclq,dclz 8-11 

dlyap2 8-12 

dnwthom 8-13 

dqcosts 8-14 

dstable 8-15 

eigpq 8-16 

nwtpar 8-17 

nwtprint 8-17 

rnormal 8-23 

8-: 26 

valprint 8-27 
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8.1 Commands Grouped by Function 


Initialization Routines 

balcred 

balanced controller reduction 

diqg 

discrete Unear quadratic gaussian design 


Homotopy Algorithm 

dnwthom 

discrete Newton homotopy algorithm 

ntwpar 

set default parameters for dnwthom 

nwtprint 

print parameters for dnwthom 

valprint 

print algorithm run-time statistics 


Controller Bases 

ccf 

convert to controllable canonical form 
(valid only for SISO controllers) 

rnormai 

convert to real normal modal form 
(a special case of tridiagonal form) 

rpf 

convert to real (or second-order) polynomial form 
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Costs 


dqcosts 


discrete costs 


Closed-Loop Matrices 

cla 

construct state matrix 

clr 

construct performance weight 

civ 

construct noise intensity or covariance 

dclp 

construct discrete observability grammian 

dclq 

construct discrete controllability grammian 

dclz 

construct discrete Z matrix 
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Utility Routines 

beam 

provide plant, noise statistics and perform weights 
for a simply-supported beam 

bodeplot 

plot magnitude and phase on same screen 

c2dv 

discretize disturbance intensity 

dlyap2 

discrete Lyapunov solution using diagonal basis 

d stable 

determine discrete stability 

eigpq 

ordered eigenvalues of product PQ 

tol80 

converts phase vector to lie in interval [- 180,180] 
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balcred— balcred 

Purpose: 

Compute a reduced-order controller using balanced controller reduction. 

Synopsis: 

[ac,bc,cc] = balcred(acO,bcO,ccO,pgram,qgram,nc) 

Description: 

Computes a controller (A c ,B c ,C c , D c ) of order n c given an initial controller of greater di- 
mension (A c fl, B C fl,C C fi, D Cj o) and the corresponding observability and controllability grammians 
i^Pgram and Q gram)- 

See also: 
dlqg 


8-6 


October 1993 


GASD-HADOC 


Harris Corporation 


beam beam 

Purpose: 

Compute a continuous- time or discrete-time representation (including noise statistics and per- 
formance weights) of a beam example. 

Synopsis: 

[a,b,c,d,rl,r2,rl2,vl,v2,vl2] = beam(nmodes,h) 

Description: 

Computes a continuous-time or discrete-time representation of the beam example described in 
the following reference: 

D.S. Bernstein, L.D. Davis and D.C. Hyland, “The Optimal Projection Equations for Reduced- 
Order, Discrete-Time Modelling, Estimation and Control,” AIAA J. Guid. Contr. Dyn., Vol. 
9, pp. 288-293, 1986. 
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bodeplot bodeplot 

Purpose: 

Plot magnitude and phase information on same screen. 

Synopsis: 

bodeplot (fhz ,m ag ,ph ase ,ti tlen ame , axes ) 

Description: 

Plots magnitude and phase on subplots that appear on the same screen. If axes is present, it 
is the 2x4 matrix of axis limits for the magnitude and phase (i.e., axes = [axismag; axisphase]. 
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ccfL _ ccf 

Purpose: 

Transform a single-input, multi-output system to controllable canonical form. 

Synopsis: 

[a,b,c,T,Tinv] = cc^aO^OjCO) 

Description: 

Transforms a single-input, multi-output plant to controllable canonical form and also returns 
the transformation matrix and its inverse. 
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cla,clr,clv cla,clr,clv 

Purpose: 

Construct closed-loop matrices. 

Synopsis: 

acl = cla(a,b,c,d,ac,bc,cc,dc) 
rcl = clr(rl,r2,rl2,c,cc,dc) 
vcl = clv(vl,v2,vl2,b,bc,dc) 

Description: 

Function cla computes the closed-loop state matrix using (2.7). Function clr computes the 
closed-loop performance weight using (2.9)-(2.10). Function civ computes the closed-loop distur- 
bance intensity or covariance using (2.14)-(2.15). 

See also: 
dclp, dclq 
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c2dv____ c2dv 

Purpose: 

Discretize a continuous- time disturbance intensity matrix. 

Synopsis: 

vd = c2dv(v,a,h) 

Description: 

Converts a continuous- time disturbance intensity V into an equivalent discrete-time covariance 
(assuming a zero-order hold with sample period h) using 


Vd 


-s: 


exp(A s )V exp( A J )ds. 
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dclp, dclq, dclz dclp, dclq, dclz 

Purpose: 

Compute the discrete, closed-loop grammians and Z matrix. 

Synopsis: 

pci = dclp(acl,rcl) 
qcl = dclq(acl,vcl) 
zcl = dclz(pcl,qcl,acl) 

Description: 

Function dclp returns the closed-loop discrete observability grainmian satisfying 

P c t = A^fPclAd + R c l- 

Function dclq returns the closed-loop discrete controllability grammian satisfying 

Qcl = AciQctAji + Vd. 

Function dclz requires the outputs of dclp and dclq to return the closed-loop discrete Z matrix 
satisfying 

Zcl = QclAjfPct- 

See also: 

cla, clr, civ 
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dlyap2 dlyap2 

Purpose: 

Solve the discrete-time Lyapunov equation by transforming to the modal basis. 

Synopsis: 

q = dlyap2(a,v) 

Description: 

Computes the solution Q to the discrete-time Lyapunov equation 

Q = AQA 1 + V 

by transforming to the complex modal basis. If the input A is a column vector, then the system is 
assumed to be in the diagonal basis and the eigenvalues are the elements of A. 
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dnwthom dnwthom 

Purpose: 

Compute am optimal discrete-time controller using the Newton homotopy algorithm. 
Synopsis: 

[ac,bc,cc,dc,par,val] = dnwthom(acO,bcO,ccO,dcO, . . . 

a0,b0,c0,d0,rl0,r20,vl0,v20,vl2, . . . 
af ,bf ,cf ,df ,r 1 r ,r2f , v If ,v2f,v 1 2 ,par) 


Description: 

Computes an optimal discrete-time controller using the Newton homotopy algorithm described 
in Section 6. On input, the vector par contains the variable algorithm parameters whose default 
values are set using function nwtpar as follows: 

par = nwtpar. 

See the nwtpar reference pages for a detailed description of the elements of par. On output, val 
is a vector containing descriptions of important run-time parameters. In particular, 


val(l) = 
val(2) = 
val(3) = 
val(4) = 
val(5) = 
val(6) = 
val(7) = 
val(8) = 

See also: 


value of homotopy parameter on return 
total number of Hessian calculations 
min # of corrections for fixed lambda 
max # of corrections for fixed lambda 
minimum homotopy step size 
maximum homotopy step size, 
number of mega-flops required for run 
number of seconds required for run. 


nwtpar, nwtprint 
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dqcosts_ _ — — dqcosts 

Purpose: 

Compute each of the quadratic costs for a given discrete-time system. 

Synopsis: 

[costs, p22, q22] = dqcosts(ac,bc,cc,dc,a,b,c,d,rl,r2,rl2,vl,v2,vl2) 

Description: 

Computes the quadratic costs for the given discrete-time system. On return costs is a 5th order 
vector whose elements have the following values: 

costs(l) = state cost ( x t Rix ) 

costs(2) = input cost (u T R 2 U) 

costs(3) = cross cost (2x T Rnu) 

costs(4) = feedthrough cost (tr Dj R^D^) 

costs(5) = total cost (sum of the above). 

The matrices p22 and q22 are respectively equal to the (2,2) blocks of the closed-loop observ- 
ability grammian ( P c i ) and controllability grammian ( Q c e )* If the controller is an LQG controller, 
then p22 = P and q22 = Q. 

See also: 

dpcl, dqcl, dcost 
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dstable — dstable 

Purpose: 

Determine the discrete-time stability of a matrix. 

Synopsis: 

sflag = dstable(a) 

Description: 

Determines the stability of the matrix A in the discrete-time sense (i.e., are the eigenvalues 
of the matrix in the closed unit circle). On return, sflag = 1 if A is stable and sflag=0 if A is 
unstable. 
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eigpq eigpq 

Purpose: 

Return the ordered eigenvalues of the product of two input matrices P and Q. 

Synopsis: 

eigpq(P,<?) 

Description: 

Computes and prints the ordered eigenvalues of the product of two input matrices P and Q. 
If the matrices are controller grammians (i.e., P and Q) the ordering can be used to determine the 
order of a reduced-order controller. 

See also: 

balcred 
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nwtpar nwtpar 

Purpose: 

Set the default parameters for the Newton homotopy algorithms. 

Synopsis: 

par = nwtpar 

Description: 

Sets the variable algorithm parameters for the homotopy algorithms for optimal, discrete-time, 
reduced-order controller design. A description of each of these parameters is given in the following 
table. 

See also: 

nwtprint, dnwthom 
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ALGORITHM PARAMETERS 


No. 

Function 

Default 

Description 

1 

Print Option 

i 

Controls amount of input during the 
execution of the homotopy algorithm 

2 

Prediction 

Tolerance 

1.0e+3 

Maximum allowable gradient norm for 
prediction step. 

3 

Correction 

Tolerance 

1.0e-4 

Maximum allowable gradient norm for 
intermediate correction loops. 

4 

Final 

Tolerance 

1.0e-4 

Maximum allowable gradient norm for 
the final correction loops. 

5 

Minimum 
Step Size 

1.0e-6 

Minimum allowable step size of the 
homotopy parameter. 

6 

Maximum 

Corrections 

10 

Maximum number of correction loops 
for a fixed value of the homotopy. 
parameter 

7 

Hessian for 
Prediction 

0 

0 uses Hessian from last correction step 
for prediction. 1 computes a new Hessian, 
parameter 

8 

Line Search 

0 

0 does not use line search 
unless cost is not decreased. 

1 always uses line search. 

9 

Automatic Run 

1 

0 lets program run interactively. 

1 lets program run automatically 


Step Size 

.01 

On input, initial step size 
On output, last step size used 

11 

Controller Basis 

1 

1-no basis 
2-tridiagonal form 

3- second-order polynomial form 

4- controllability canonical form 

12 

Prediction Option 

1 

0-no prediction 
1-linear prediction 

2- circular arc prediction 

3- cubic spline prediction 
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nwtprint 

Purpose: 

Print the variable homotopy parameters. 

Synopsis: 

par = nwtprint(par) 

Description: 

Prints the information contained in the vector par that describes 
rameters for the Newton homotopy algorithms. 

See also: 

nwtpar, dnwthom 
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rnormal rnormal 

Purpose: 

Convert a plant to real normal modal form or a standard alternative. These forms are special 
cases of the tridiagonal form. 

Synopsis: 

[at,bt,ct,dt] = rnormal(a,b,c,d) 
or 

[at,bt,ct,dt] = rnormal(a,b,c,d, ‘modal’) 

Description: 

The first call with four input arguments converts a plant to real normal modal form, i.e., the 
transform of state matrix, at has 2x2 blocks of the form 

-V{ Wj 
-Vi 

If ‘modal’ is input as a fifth input argument, that at has 2x2 blocks of the form 

0 1 * 

-(<T 2 + LJ 2 ) 2(7 

See also: 

trimats 
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rpfx rp fx 

Purpose: 

Transform a system to second-order polynomial form. 

Synopsis: 

[A,B)C,T,Tinv] = rpfx(A0,B0>C0) 

Description: 

Transforms a system to second-order polynomial form and also returns the transformation 
matrix and its inverse. 
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tol80 to 180 

Purpose: 

Transform a phase vector to lie in the interval [-180,180]. 

Synopsis: 

phasel80 = tol80(phase) 

Description: 

Transforms a phase vector to lie in the interval [-180,180]. 
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valp ri nt valprint 

Purpose: 

Print the run-time homotopy parameters. 

Synopsis: 

valprint(val) 

Description: 

Prints the information contained in the vector val that describes important run-time informa- 
tion for the Newton homotopy algorithms. 

See also: 

dnwthom 
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Appendix A: Cost Derivatives 

In this appendix we consider the cost function J(A c ,B C yC c ) defined by (2.4) or, equivalently, 
(2.16). We derive expressions for and 

Let C denote the Lagrangian defined by (2.33) which is rewritten here as 


£(A C , B cy C c * 

Z? c ,P,Q) = tr[gP + P(A<?A T 

+ V-Q) + DjR 2 D c V 2 ). 

(A 2) 

where 

Q = AQA t + V 


(A.3) 

Then, 




II 

V 05 

dj dj dj dj 

dB c ~ dB c ' dC c ~ dC c ' 

^ cS 

II 

. y 

b c* 

(A.4a,6,c, d) 

subject to the constraint 

dQ 


(A.5) 

or, equivalently, 

P = A 1 PA + R. 


(A.6) 


Now, let <t> denote an element of A c ,B c ,C c or D c . Then, 


3J_ _dC 

d<f> ~ d<{> 


or equivalently 


where 


= tr 


- 9R T :^9A T dV x 

® d<j> + P ^d(f>^ A +A ® d<t> + 


+ ^tr (D c rt R 2 D c V 2 ) 


% - "Wf + p ^ + 2 ^ ) + ^ tr(B * 


Z = Q/TP. 

It follows from (A.8) that 

dj dK dj dK dj dK dj dK 
dA c ~ dA' dB~dB' dC c ~dC' dD c ~dD c 


where 


K(A e , B c , C c , D C ,P,Q, Z) ^tr [QR(C c , D c ) + PV(B C , D c ) 

+2A(A c ,B c ,C c , D C )Z } + tr (Dj R 2 D c V 2 ) 


(A.7) 

(A.8) 

(A.9) 

(A. 10) 

M-ii) 
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and from (2.7), (2.9)-(2.10), and (2.14)-(2.15) 


A(B c ,C c , D c ) 
R(B C , C c ,D c ) 


A - BD c C -BC c 
B c C A c — B c DC c \ ' 

Rw R\i 

R-ll R-22 . 


where 


Rn =Ri~ C T DjRj 2 - R u D c C + C T DjR 2 D c C 
R 22 — —R\ 2 C c + c t dJr 2 c c 
R 22 = cJr 2 c c , 


(A. 12) 
(A.13) 


(A. 14a) 
(A. 146) 
(A. 14c) 


and 


where 



(A. 15) 

V n = Vi - BD C V& - V l2 DjB T + BD c V 2 DjB J 

(A. 16a) 

V 12 = VnBj - BD c V 2 Bj 

(A. 166) 

v 22 = B c v 2 Bj. 

(A. 16c) 


The desired derivative expressions will be derived using (A. 10). The development proceeds by 
considering each of the four terms in the right hand side of (A.l) and differentiating these terms 
with respect to A c ,B c ,C c and D c . It is assumed that P,Q and Z are partitioned conformably with 
A,R and V such that 


P = 


P\\ P\ 2 

,A t 2 P22 


Q 


“ [q t 


11 


12 


Ql2 

Q22 J 


Z\\ Z\ 2 

Z 2 \ Z 22 j 


(A.17) 


tr QR 


tr QR = tr(<3ii^n + 2Qi 2 R,J 2 + Q 22 R 22 ). (A. 18) 

Using (A. 14), we may then write 

tr QR = tr(g n fli - 2 Q 11 C t DJrJ 2 + Q u C T DjR 2 D c C ) 

+2tr(-Q 12 C T ^ 2 + Q 12 CJR 2 D c C) (A.19) 

+ti(Q 22 Cj R 2 C c ). 
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Differentiating (A. 19) gives 

sb* 

^-UQR = 0 

■^-tiQR = 2(—Rj 2 Qi2 + RjDcCQu + RiC c Q 22 ) 

^tr QR = 2(-Rj 2 Q u C T + R 2 D c CQ n C T + R 2 C c Qf 2 C T ). 

tr PV 


tr PV = tr(AiVii + 2 Aa<S + A2V22). 

Using (A. 16), we may then write 

tr PV = tr(P u Vi - 2P n BD c V& + P u BD c V 2 Dj B t ) 
+2tT(P u B c V? 2 -P 12 B c V 2 D]B T ) 

+U(P 22 B c V 2 Bj). 

Differentiating (A. 25) gives 

d 


dA c 

d 

dB c 

d 

dC c 

d 

dD c 


tr PV = 0 

trPV = 2{P? 2 V U - P? 2 BD c V 2 + P 22 B c V 2 ) 
tr PV = 0 

tr PV = 2(-B T P u V 12 + B t P u BD c V 2 - B r P 12 B c V 2 ). 
tr AZ 


tr AZ = tr[(A - BD c C)Z n - BC c Z 2l + B c CZ 12 + (A e - B c DC c )Z 22 \ 
Differentiating (A. 30) gives 

d - - ~ T 

— - trAZ = 2 Z?>. 

V Ac 

T^j-tTAZ = 2(Zj 2 C T - Zj 2 CjD T ) 

^tAZ = -2(B T Zjy + D T BjzJ 2 ) 

- d ^-UAZ = -2B T Zj 1 C T . 
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(A.22) 

(A.23) 

(A.24) 

(A.25) 

(A.26) 

(A.27) 

(A.28) 

(A.29) 

(A.30) 

(A.31) 
(A.32) 
(A. 33) 
(A- 34) 
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tr DjRjDc}^ 


R 2 D c V 2 = 0 

(A.35) 

JL tT DjR 2 D e V 2 = 0 

(A.36) 

^tr DjR 2 D c V 2 = 0 

(A.37) 

J-XtD1R 2 D c V 2 = 2R 2 D c V 2 . 

oD c 

(A.38) 


«/(A ct i? c ,C c .D e ) 

It follows from (A. 10) and (A. 11) with (A.20)-(A.23), (A.25)-(A.29), (A.31)-(A.34) and 
(A.35)-(A.38) that 

f£ ■ « < A - 39 > 

= 2(P 1 T 2^12 - PuBD c V 2 + P 22 B c V 2 

+ Zj 2 C T - Zj 2 Cj D t ) (A.40) 

|£ = 2(-Rj 2 Qu + R 2 D c CQ 12 + R 2 C c Q 22 
oC c 

- B T Zj 2 - D r BjzJ 2 ) (A.41) 

= 2(-RJ 2 Q u C t + R 2 D c CQu + R 2 C c Qj 2 C T 

o D c 

+ B T P„ v 12 + b t p u bd c v 2 - B T P l2 B c V 2 

- B t Z^C* 1 + R 2 D c V 2 ). (A.42) 
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Appendix B: Closed- Loop Matrix Derivatives 


In this appendix we present explicit expressions for the derivatives 


where 


9 = 


" vec(i4 c )' 
vec(5 c ) 
vec(C c ) 
,vec(D c ). 


- A- BD c C -BC c 
4_ B c C A c - B c DC c ’ 


R = 


it n 

rJ 2 


it 12 

it 22 


where 


R n = Ri - C T DjRj 2 - R u D c C + C T DjR 2 D c C 

Ri2 = ~R\iC c + c t dJr 2 c c 

R-22 = C]R 2 Cc, 


and 



V12 

V 22 


where 


Vii = V 2 - BD c V? 2 - V u DjB r + BD c V 2 DjB T 
Vn = Vi 2 B] - BD c V 2 Bj 


v 22 = Bjv 2 Bj. 


It is assumed that the plant matrices (A,B>C, D) y the cost weighting matrices (i?i, 
and the disturbance matrices (Vi, V^V^) are the following functions of A. 


where 


A(X) 

C( A) 


B( A) 
D(A)J 


A o 
Co 



Bj 

Df J 


Aq Bo \ 
Co Dq _ ) 


RU A, t(A)] “ 


Lr(X) = Lfifl + A (Lrj - LRfi) 


(B.l) 

(B-2) 

(B.3) 

(B.4a) 

(B.46) 

(B.4c) 

(B.5) 

(B.6a) 

(B.6h) 

(B.6c) 

^12 >^2) 

(B.7) 

(B.8a) 

(B.86) 
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and Lfi t o and Lrj satisfy 


where 


L «,oiio = [$,•; 

LrjLr f = \ jp ,f 


Ru,o 

Rifi ' 

R\2J 

Rij . 


Vi (A) Vn(XY 

LKS(A) ^ T (A)J 


= LvWLjr(X) 


Lv( A) = Lv,o + X(Lvj — Lv, o) 


and Ly , o and Xvj satisfy 


- r T 

y Vo^V t O 


rT 



Xo 

W 


. ViT.o 

V 2 , 0 . 



Vi 2./' 


. ^12,/ 

Vi.0 . 


Below, we use the notation 


Note that from (B.7)-(B.9) 




dX 


■ A B 


Aj — Ao Bj — Bq 

C D 


Cj — Co Dj — Dq 


where 


where 


Ri R\2 

R-'i : 2 R-2 


= L r L r 


Lr — Lrj — Lra 


Vi V X2 

r 
2 


Vr T 2 V, 


= Lv Lv 


Lv = Lvj — Lv, o- 


(B.8c) 

(B.8d) 

(B.9a) 

(B.96) 

(B.9c) 

(B.9d) 

(B. 10) 

(Bll) 

(B.12a) 

(B.126) 

(B.12c) 

(B.12d) 


The derivations of the expression for jf- , and are primarily based on the application 
of the following derivative formulas. It is assumed that X is an n X m matrix. 
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where 


Derivative Formulas 


-^-XA = [A(j,:)U^-i 

axij 

(B.14) 

~r~~AX = [A(:,t)] co |_j 

dx{j 

(B.15) 

4- ^ A =l A (i, Olrow-i 

dx{j 

(B.16) 

-^-AX T = [A(:,j)U-i 

(IX ij 

(B.17) 

■y—AXB = A(:,i)B(j,:) 
ax j j 

(B.18) 

-/-AX r B = A(:J)B(i,:) 

dx.j 

(B.19) 


dA 

dfj 


dA 

do-c,kt 


dA 

db Ct kt 

dA 

dCc^kl 


0 0 

o !«&.. 

0 0 

[C(^i : )]row — it [DCc (£, ')lrow-fc 

0 -[B(:,fc)] co i_* 

0 #*--[B e D{:,k)U-t 

aCc k,t J 


-B(:,k)C(i,:) 0 

0 0 

Q b A * t and are given respectively by (D.2.36) and (D.2.37) of Appendix D. 

dR 

de. 


dA 

dd Ci ki 


(B.20) 

(B.21) 


(B.22) 


(B.23) 


dR 

da c< k( 

dR 

db C 'kt 


(B.24) 

(B.25) 
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o 


[-An(:,fc) + C T i>Ji?2(:,fe)]coi-< 


dk 

dc c ,ki 


dR 

dd c ,ki 


SYM [CjR 2 (:, k)U- t + l&Cc jk, :)] row _, 


‘ C T (:,l)[R2DcC(k,:)-Rj 2 (k,:)] 
^CjDjR 1 (:,k)-R l 7 (:,k))C(t,-.) 

SYM 


C T (:,l)fl 2 C c (fc,:) 


0 


dV 

dOj 


dV 

d<icM 


= 0 


dV 

db cM 

dV_ 

dc C 'ke 

dV 

dd^kt 


0 [ V n (:,e)- BDcV 1 (:J)U. k 

SYM [F 2 5j(^:)] row _ fc + [gcjjQ.^cd-fc 


= 0 


B(:.fc)l V 2 DT B T (/;) _ K T (<:)| 

+[£DcV2(:,/)-V, 2 (:/)]flT(fc i: ) 

SYM 


0 




dA 

M. 


- _ [ A - BD c C - BD c C 
B c C 


-BC c 
— B c DC c 


R± 


dR 
d A 


where 


R = 


iZn R\2 
:,T :T 

R\2 ^22 


k u = Ri - C r DjRl j - C t DJR? 2 - R l2 D c C - R l2 D c C 
+ C T DjR 2 D c C + C T DjR 2 D c C + C t dJr 2 D c C 
k n = -R u C c + C T DjR 2 C c + C T DjR 2 C c 
r 22 — C]R 2 C c . 


(B.26) 

(B.27) 


(B.28) 

(B.29) 

(B.30) 

(B.31) 


(B.32) 


(B.33) 


(B.34a) 

(B.346) 

(B.34c) 
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where 



'^n 

fn ^ 22 . 


= Vi - BD c V ? 2 - BDcVyi - VuDjB 7 - V l2 D? B t 
+ BD c V 2 D c B t + BD c V 2 D c B t + BD c V 2 D c B r 
Pu = -V n Bj - BDjV 2 Bj - BD c V 2 Bj 

v 22 = b c v 2 bJ. 


(B.35) 


(B.36o) 

(B.366) 

(B.36c) 
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Appendix C: The Input-Normal Riccati Basis 

The homotopy function H{9 , A) described in Section 4 is defined to solve the optimal reduced- 
order dynamic compensation problem for discrete-time systems. The elements of the vector 0 
include parameters which completely describe the controller (A c , R c ,C c ,I} c ). For computational 
efficiency it is desired that the vector 0 be as small as possible. Thus, we desire to represent the 
controller (A c , I? c , C c , D c ), assumed to be minimal, with the fewest parameters possible. The results 
of this section reveal that in a certain basis, which we will denote as the input normal Riccati basis, 
the controller plant matrix A c is almost always completely characterized by its input and output 
matrices B c and C c . 

Theorem C.l. For every minimal n* h order plant (A c ,B cy C c ,D c ) there exists a similar- 
ity transformation Tj and a positive matrix S7 = dia.g{c*^ t } Jtlj such that ( A c = Tj~ l A c TjB c = 
TJ 1 B C , C c = C c Tj, D c — D c ) satisfies 

0 = A c + A? + B c Bj-C]C c (C.l) 

0 = A*n + nA c + cJc c -nB c Bjn. (c.2) 


In addition, 


= [(CjCch/iBcBjU '' 2 
A cM = \[(CjC c ) ii -(B c Bj) tt ] 

and if a;,- ^ a )j for i ^ j, then 

A • • — ^j(l + ~ (CjC c )ij( 1 + U?j) ^ ^ j 

CyiJ U?i — Uj 

so that A c is completely and uniquely determined by B c and C c . 


(C.3) 

(C.4) 


(C.5) 


Proof. The minimality of (A c , B c ,C c ) insures that there exists unique positive definite solu- 
tions Q and P of 


0 = A C Q + QAj + B c Bj - QCjCcQ (C.6) 

0 = AjP + PA c + Cj C c - PB c Bj P. (C.7) 

It is well known that there exist a transformation Tj such that 

T^QTj r = I nc , TjPTi = (l (C.8) 
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It follows from (C.6)-(C.8) that (A c ,B c ,C c ) satisfies (C.l) and (C.2). 

We now show by construction that (C.3)-(C.5) hold. First recognize that (C.l) and (C.2) are 


equivalent respectively to 

0 = Ac, a + Ac,a + {BcBDa - (C c T C c ) 0 (C.9) 

0 = UiA c ,ij + u>jA c ,ji + {CjC c ) tJ — u)iUj(B c Bj)ij. (C.10) 

Letting * = j in (C.6) gives (C.4) while it follows from (C.7) that for i = j 

0 = 2uiA c ,a + (Cj C c )u — <*>i(Bj B c )u. (C.ll) 

Substituting (C.4) into (C.ll) gives 

0 = -(BjBcU^ + [(CjC c )ii - (BcBj^Pi + (CjC c )u (C.12) 

which has positive solution given by (C.3). 

Multiplying (C.9) by -u- and adding the result to (C.10) gives 

0 = (u>i — u>j)A C 'ij + (CjC c )y(l + uj) - Uj{ 1 + Ui){B c Bl ) t j (C.13) 

which implies that if u>i ^ Uj for i ^ j then A Ct ij is given by (C.5). □ 


Definition C.l. If the minimal order plant ( A c ,B cy C c ,D c ) satisfies (C.l) and (C.2) of Theo- 
rem C.l, then the plant is said to be in input normal Riccati form . 

Remark C.l. Input normal Riccati form is similar to the input normal form of Moore [1] 
which is further explored by Kabamba [2]. 


Now, define 

ft = diag{u>j}".f j 

(C.14) 

and H € such that 

h .. ± 1 ( w < = 

l o, t = j. 

(C.15) 

or equivalently 

ft ^ diag(C c T C c )diag(5 c 5j)- 1 

(C.16) 


H = (N nc - /„J/[iV nc ft - SlN nc + J n J 

(C.17) 
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where » m € IR mxm defined by 

N miij = 1. (C.18) 

Then the following remark holds. 

Remark C.2. If Uj for i j, then (C.4) and (C.5) are equivalent to 

A e = ±[CjC c - B c Bj] * I nc + [CjC c (I + ft) - ( I nc + Sl)B c BjSl ) * H (C.19) 

Proposition C.l. Let A and Z be in IR" X " with A diagonal. Then, 

AZ = A*Z (C.20) 

and 

ZA = Z*/F. (C.21) 

Remark C.3. It follows from Remark C.2 and Proposition C.l that if u>, ^ u>j for i ^ j, then 
A c can be computed by 

Ac = l[CjC c -B c B])*I nc +[CjC c +(C c Cj)*n T -(B c Bj+Sl*(B c Bj+n*(B c Bj))*n r ]+H (C.22) 

where from (C.17) with (C.20) and (C.21) 

H = (N nc - I nc )/fi T - ft + I Jlc ). (C.23) 
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Appendix D: The Gradient of the Cost Functional for the Input Normal Riccati Basis 

In this Appendix it is assumed that the controller (A c , B c ,C c ,D c ) is in the input normal 
Riccati form of Appendix C and is hence completely described in terms of B c> C c and D c . We 
let j(B c ,C c ,D c ) be the restriction of the cost functional J(A c ,B c ,C c ,D c ) defined by (2.4) or 
equivalently (2.16) on the set of generic input-normal Riccati triples (A c , B c ,C c ). Also, define 

vec(R c ) 

9 = vec(C c ) . (D.l) 

vec(D c ) 

Then, with some abuse of notation we can write the restricted cost functional as J{0). The 
homotopy algorithms to be defined later will be based on finding 6 satisfying 

° = /(*) = (D2) 


Now, recognize that 


vj(«) T i °4 = 


(D.3) 


The next theorem present very useful expressions for and This result is very similar 

to Theorem 2 of [D.l]. 

Theorem D.l. The derivatives and are given by 


wrti +2(Y - nzn)B ‘ 

^ = g + 2C t( z-y) 

dj dj 
dD c dD c 

where Y 6 IR nxn and Z G IR nXn are symmetric and satisfy 

0 = |^ + 2(7 + QZ) 

U Aq 

0 = {{A c - BcBltyZ)^ i = l,2,...,n c . 


(D.4) 

(D.5) 


(D.6) 


(D.7) 


Proof. Since the triple (A c ,i? c ,C c ) is constrained to be in input-normal Riccati form, it 
satisfies 



0 = A c + A t c + B c Bj - CjC c 

(D.9) 


0 = Ajn + ClA c + CjC c - SlB c Bjn. 

(D.10) 
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Following the proof of Theorem 2 of [D.l] we define the new Lagrangian 

£(A c ,B c ,C c ,D c ,P,Q,Y,Z) 

= C(A c ,B c ,C c ,D c ,P,Q) (DU) 

+ tr \Y(A e + A?- CjC c + B c Bj) + Z{A T C S 7 + S2A C + CjC c - SIB c B]Q)} 


where Y and Z are n x n symmetric Lagrange multiplier matrices. Then 



dj 

dj dj dj dj dj 

(D.12a,6, c) 


dB c 

~ dBj dC c ~ dCj dD c dD c 

subject to the constraints 





„ 9C n dC 
dAj dn' 

(D.13a,6) 

Now, 


ei = H +2(Y + nz) 

(D.14) 

and 

dC 

3Q 

= A C Z + ZA t c - ZQB c Bj - B c BjQZ 

(D.15) 

which implies 





dC 

dui{ 

— 2(A C Z — B c Bj fiZ)u , i — 1, . • . , n c . 

(D16) 


Equations (D.7) and (D.8) follow respectively from (D.13a) with (D.14) and (D.13b) with (D.15). 
Finally, (D.4)-(D.6) follow respectively from (D.12a,b,c). □ 


We now state a very important corollary which describes how to efficiently compute Y and Z 
satisfying (D.7) and (D.8). For convenience, we define 

i (D.ira) 

F=A c -B c B t c Q. (D.175) 


and rewrite (C.16) here as 


D i diag(C c T C c )diag(5 c B c T )- 1 . 


Corollary D.l. The matrices Y and Z in (D.7) and (D.8) satisfy 

Y = -(\c A . +nz) (D.19) 

i = i 
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and 

o = zn-nz + ±{c A <-c T Ac ) (D.22) 

which if u>i / Uj for » / j is equivalent to 

= (D.23) 

Uj — 

Proof. Equations (D.19) and (D.20) follow respectively from (D.7) and (D.8). Since Y is 
symmetric 

0 = y-y T . (D.24) 

Substituting (D.19) into (D.24) gives (D.22) or equivalently 

0 = (uj - u>i)zij - ^(C Ac ii - C Ac) , ), i ? j (D.25) 

which if u>i ^ Uj for i ^ j is equivalent to (D.23). □ 

Remark D.l. If Ui ^ uj for i ^ j (D.20) and (D.21) are equivalent to 

Z = Z 0 + (D.26) 

where 

(D.27) 

Z d = -diag(F)" 1 *(FZo) (D.28) 

and H is given by (C.23), rewritten here as 

H = ( N nc - 7 n J/[ft T -ft + 7„J. (D.29) 

Expressions for the partial derivatives ■§£-, licl' an ^ are derived * n Appendix A. Here, 
we cite only the final results. First, we define 

Z = QA t P (D.30) 


and note that P, Q and Z have the partitioned forms 

Z11 I12] 

Z\2 ^22 J 


P = 


P 11 Pl2 
Fu ftd 


. 0 = 


Q 11 Qn 

1$12 O22 


Z = 


(D.31) 
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where the (1,1) and (2.2) blocks of each matrix are respectively n x n and n c x n c . With this in 
mind, the cost derivatives are given by 

f£ = «« < D - 32 > 


= 2(/j2 V 12 - PnBD c V 2 + P 22 B C V 2 
+ Z? 2 C r - Zj 2 CjD T ) 

-Jtt — 2 (~ Ri 2 Qi2 + RiDcCQu + R 2 C c Q 22 

(s(s c 

+ b t zJ 1 -d t bJzJ 2 ) 

jfi- = 2(- j R? 2 Q 11 C t + R 2 D c CQ n C T + R 2 C c Qj 2 C T 

- b t p u v 12 + b t p u bd c v 2 - b t p 12 b c v 2 


- B 1 Z 1 T 1 C T + R 2 D c V 2 ). 


(D.33) 

(D.34) 


(D.35) 
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Appendix E: The Homotopy Map and It’s Jacobian for the Normal Riccati Basis 


As stated in the previous Appendix the objective is to find 8 satisfying 

m = o (E.i) 

where 

f(6) = VJ(0) T (E.2) 

and J(0) denotes the restricted cost functional for the input-normal Riccati basis. In this section we 
define a homotopy map to accomplish this task and show how to efficiently compute it’s Jacobian. 

Definition of the homotopy map h(0, A) 


To define the homotopy map we assume that the plant matrices (j4, 5, C, D), the cost weighting 
matrix (i?i, # 2 ^ 12 ) and the disturbance matrices (V \ , V 2 > ^ 12 ) are functions of the Homotopy 
parameter A £ [0, 1]. In particular, it is assumed that 
A( A) B(\y 


[ RiW 

r ^(A) 

[v®A) 


C( A) 0(A) J 


Aq Bo 

Co Do 


+ 


* 12 (A) 
* 2 ( A) 

Vj 2 (A) 

W) J 


*1,0 

pT 

^12,0 


*12,0 
*2,0 

Vi l0 ^12,0 

1 * 12.0 ^ 2,0 


+ A( 


(\A } 

*/]_ 

Aq Bq 

h 


*>/J 

Co Do 

V 

* 1 ,/ 

* 12 ./] 

_ * 1,0 

* 12,0 ] 

DT 

." 12 ,/ 

* 2 ,/ J 

.* 12,0 

* 2,0 J 

V 1J 

^!2./l 

_ ' F ,, 0 

Fl2 , o ]\ 

* 12 ,/ 

^ 2 ./ J ' 

.* 12,0 

F 2 ,o JJ 


(E-3) 

(E-4) 

(E-5) 


Note that (E.3)-(E.5) imply that .4(0) = A 0 and A( 1) = Aj, B( 0) = B 0 and 5(1) = Bf, etc ... 
and it is understood that Aj,Bj , . .. were referred to in the previous sections simply as A,B r . ... 
The change in notation is simply for convenience. 


The homotopy map h(0, A) is defined by 


where 


MM) 


-vec (H Bc (8,X))' 
vec (5 Ce (M)) 


HbS8, A) = 2(5^12 - P? 2 BD C V 2 + P 22 B c y 2 

+ z 1 t 2 c t - zJ 2 cJd t + (f - nzsi)B c ) 
Rc c {8, A) = 2(— RJ 2 Qu + RiD c CQn + * 2 C c Q 22 

- P T Z£ - d t bJzJ 2 + C c (Z - F)) 
Rd c (8, A) = 2(-P? 2 0iiC ,t + R 2 D c CQ n C T + 

- 5 t AiF ]2 + B r P n BD c V 2 - B t P 12 B c V 2 

- B t Z j.C' 1 + R 2 D c V 2 ). 


(E.6) 


(E.7) 


(E.8) 


(E.9) 
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Here, it is assumed that P,Q and Z satisfy 

P = A r PA + R 

(E.10) 

Q = AQA t + V. 

(E.11) 

Z = QA t P. 

(E.12) 

In addition, Y and Z are given by 

Y = + (IZ) 

(E.13) 

z = z 0 + z d 

(E.14) 

where 

Zo = \(Ca.-C t Ac )*H 

(E.15) 

Z d = - diag(F)- 1 * (FZo) 

(E.16) 

C Ac = QZj 2 

(E.17) 

F = A c - B c B]Sl 

(E.18) 

fl ± diag(CjC c ) diag(5 c 5j) -1 

(E.19) 

H = {N nc -I nc )./[Q. T -n + 1 

(E.20) 

Note that (E.13) and (E.14) are equivalent to 

Q = £ Ac +2(Y + QZ) 

(E.21) 

0 = [F Z]ii , t = 1,2, .. . ,n c . 

(E.22) 

Also, note that it follows from the results of the previous section that 

MM) = /(*)(= vA*) T )- 

(E.23) 


Also, note that h(0, A) is the transposed gradient for the optimization problem with parameters 
We now consider that computation of th e Jacobian of h(0,\). Note that 


, ,. T \dh dll'] 

= [as ^aJ- 


(E.24) 
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Recalling that 9 is defined by (D.l), such that for some integers k and l, 9j is given by 


@j — bc,k ( » 9j — c Ct ic t , or 9j — ^c,ki • 


(E.25) 


It follows from (E.6) that is of the form 


dh 

d9 


■ • •' vec (efe^) • • * ve <dij; H Bc ) ■ ■ ■ vec( ) 
• • • vec ( ) • • * vec ( d£j ; Hc < ) • • • ve <d£j ; Hc ‘ ) 

9 


L • • • vec (e 537 ^J • • •’ vec(-s£-H Dc ) . . . vec( e^H Dc ) J 


(E.26) 


and ^ can be expressed as 


dh 

d\ 


vec (jx H Bc) 


i T 


vec(*»c.) | • (E.27) 

Below, we develop explicit expression for the derivative terms appearing on the right hand sides of 
(E.26) and (E.27). We use the notation 


M i 


DM 

dOj 


M * 9M 

M = W 


(E.28) 


Differentiating (E.10)-(E.12) with respect to 9j yields 


P U) = A T P U) A + (A U)T pa + A T PA (j) + r U) ) 

(E.29) 

qW = AQ (j) A T + (A U) qA t + AqA u)T + v u) ) 

(E.30) 

z (J) = q U) A t P + qA U)T P + qA t P u) 

(E-31) 


where expressions for the derivatives A^\R^^ and are given by (F.20)-(F.28) of Appendix F. 
Differentiating (E.21), (E.22), (E.17) and (E.18) with respect to 9j yields 


0 = (rJJ + 2 + n (j) Z] + 2 (Y U) + DZ (j) ) (E.32) 

-[F (i) Z]u = [FZ^)n, i = 1,2, . . . ,n c (E.33) 

where 

Cf, = QZif (E.34) 

F U) = A U) _ (B c Bj) u) Sl - B c Bjn {}) (E.35) 
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and the derivatives Ac*\{B e Bj)M and may be computed using the results of Appendix G. 
Note that if we define 

C = £% + 2Sl U) Z (E.36) 

then (E.32) and (E.33) are equivalent to 

yG) = + (E.37) 

Z U) - z(j) + z ( d j) (E.38) 

where 

4 J) = (E.39) 

z 

Z f = - diag(F)" 1 * (FZ ( 0 j) + F (j) Z). (E.40) 


Differentiating (E.10)-(E.12) with respect to A yields 

P = A T PA + {A PA + A 1 PA + R) (E.41) 

Q = AQA 1 + (AQA t + AQA + V) (E.42) 

Z, = QA T P + QA P + QAP. (E.43) 


where expressions for A, R and V are given by (F.29)-(F.33) of Appendix F. Differentiating (E.21), 
(E.22) and (E.17) with respect to A yields 

0 = Ca c + 2 (F + (IZ) (E.44) 

0 = lFZ]a, i = 1,2, . . . , n c (E.45) 


where 



Note that (E.44) and (E.45) are equivalent to 

Y = -(\c Ac + nz) 

Z = Z Q + Zd 


where 

Zo = \{Ca c -CI c )*H 

Zd = — diag(F) -1 * (FZ 0 ). 
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(E.47) 
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(E.49) 
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Before presenting the desired derivative expressions we define 

H' Bc (P^,Z^,Y^\Z^) = 2(P 1 ( 2 i) Vn - p[f bd c v 2 + P^B C V 2 

+ z ( 1 ? t c t - z { 2 { )T cJd t + ( y< j > - nz( j] si)B c ) (E. 51 ) 

H' C '{Q U) ,Z {i) ,Y {i \Z (i) ) = 2 (-Rj 2 Q$ + RiDcCQW + R 2 C c Q$ 

- B t Z { 2 \ )T - D T BjZ 22 T + C c (Z U) - K (J) )) (E.52) 

H' Dt (P U \Q {i \Z U) ) = 2 (-RjiQnC 1 + R 2 D c CQ^C t + R 2 C c Q[f C T ) 

- B r z[i )T C T ) (E.53) 


Notice that the right hand sides of (E.51)-(E.53) are essentially identical in form to the right 
hand sides of (E.7)-(E.9). The difference is that P,Q,Z,Y and Z and have been replaced by 
p(j) ? Q( J ), ZO’>,yU> and Z W and the last term ( c lR<iD c V'i ) in (E.9) has no counterpart in (E.53). 


Derivatives with Respect to B Cy ki 
Differentiating (E.7)-(E.9) with respect to Oj ) gives the following. 

= H' Bc (P^,Z^,Y^,Z^) + 2{P 22 E^ n V 2 + (Y - nZQ)E k( ] 

OBcM 

-2{n U) zn + nzn {i) )B c 


dHc, 

dB Ct kt 


H' c (QV\Z (j \Y( j \Z {j) ) 


2D r E { n t ' k x ) nc Zj 2 


dHi k 
dB cM 


H' Dc {P U \Q ( ’\zM) - 2B r P u E ( n k c i ) n 


V-> 


Derivatives with Respect to C c< kt 
Differentiating (E.7)-(E.9) with respect to ffj(= c Ci kt ) gives the following. 

= H' b (pU\Z ( ’\Y (i \zW) - 2ZlE ( n t ’ k) n D t - 2 + f IZQ U) )B C 

vtc'kt c c 


(E.54) 


(E.55) 

(E.56) 


(E.57) 


dH Cc 

dC Ct kt 


H' Ce (Q<»,zU\YV\zM) + 2R 2 E (k X c Qn + 2 E {k $ nc {Z - Y) 


dHp c 

dCc^t 


H' Dc (P {j) ,Q u \Z U) ) + 2R 3 E l n k J t J nc Qj t C 1 


(E.58) 

(E.59) 


Derivatives with respect to d c> kt 
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Differentiating (E.7)-(E.9) with respect to d Cyi a gives the following. 

= H' BC (PU\ZV\YM,ZM) - 2 P? 2 BE { n k y n V 2 

OVcM 

1^- = HcSQ U) » Z U) « Y (i) ,Z^) + 2R 2 E ( n k $ ny CQn 
OVcM 

1 ^- = H' Dc {P (i) ,Q {i \Z (i) ) + 2(R 2 E[ k ^ n CQ 11 C r + B T P u BEi k J l J ny 

OJJckl 


V 2 



Derivatives with Respect to A 


(E.60) 

(E.61) 


(E.62) 


Differentiating (E.7)-(E.9) with respect to A gives 

dH B , 


d\ 


= H' Bc (P,Z,Y,Z ) 

+ 2(P^V U - P? 2 BD c V 2 - P? 2 BD c V 2 + P 22 B C V 2 ) 


+ Zj 2 C r - Zj 2 Cji) T ) 

'^ = h' Cc (&2,y,z) 

+ 2(-Rj 2 Qu + R 2 D c CQn + R 2 D c CQu + R 2 C c Q 22 


dHp c 

d\ 


+ b t zJ x - d t bJzJ 2 ) 

H' Dc {pM,'z) 

+ 2(-Rj 2 QnC T - Rj 2 QnC T 

+ R 2 D c CQnC T + R 2 D c CQnC t + R 2 D c CQnC T 
+ R 2 C c CQ? 2 C t + R 2 C c Qj 2 C T ) 


+ 2(-B T P u Vn - B T PnV 12 

+ b t p u bd c v 2 + b t p u bd c v 2 + B r P 11 BD c V 2 

- B t P 12 B c V 2 - B t PuB c V 2 ) 

- 2B T Zj 1 C T - B T Z n C T 
+ 2(R 2 D c V 2 + R 2 D c V 2 ) 


where 


A B _ Aj — Aq B j — Bq 
C D Cj — Co Cf — Co 


Ri R\2 

R \2 R% 

vs V 2 J 


R\j - Ri,o Rnj - Ru,o 
.^ 12 ,/ - ^? 2,0 R-1J — ^ 2,0 _ 
v iif — Vjo Vj2,/ — Vi2,0 

y&j-vs* v 2 J -v 2fl J- 


(E.63) 


(E.64) 


(E.65) 

(E.66) 

(E.67) 

(E.68) 
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Appendix F: Closed-Loop Matrix Derivatives for the Input Normalized Riccati Basis 

In this appendix we assume that ( A c , B c ,C c , D c ) is restricted to the input-Normal Riccati basis 
and present explicit expressions for the derivatives Jgj, fjj, , fx» an ^ 7>T w here 


where 


9 = 


vec (B c ) 
vec (C c ) 
vec (D c ) 


i _\A-BD c C -BC c 

[ B c C A c -B c DC c \' 


R = 


Ru 

Rl 2 


R\2 

R-22 ’ 


(F.l) 


(F.2) 

(F.3) 


R n =Ri- C t dJrJ 2 - R u D c C + C T DjR 2 D c C (F.4a) 

R 12 = -RnC c + C T DjR 2 C c (F.44) 

R 22 =CjR 2 C c , (F.4c) 


and 

(F.5) 

where 


V = 


VJ, V u 
lV3 v„ 


V n = Vi - BD C V 3 - V n DjB r + BD c V 2 DjB T 

(F.6a) 

V 12 = V u Bj - bd c v 2 bJ 

(F.6b) 

v 22 = bJv 2 bJ. 

( F.&c) 


It is assumed that the plant matrices ( A,B,C,D ), the cost weighting matrices (Ri,R 22 ,R 2 ) 
and the disturbance matrices ( V , Vi 2 , Vi ) are the following functions of A. 


MW s(A)] 

\A 0 

Bo' 

+ a( 

Af 

B/l 


Bo' 



lc(A) B(A)J 

~ [Co 

Do. 

c f 

D/J 

.Co 

Do. 

) 


' Ri( A) RnW 


Ri y o 

Rl2,0 

+ a(| 

' Ru 

B 12 ,/ 


Bi,o 

R 12 ,< 

) 

RJ 2 (X) R 2 ( A) 


r t 

. n 12,0 

R 2 ,0 . 

R r 

l n 12J 

f? 2 ,/ . 


. B?2,o 

R 2 ,0 

- 

■ Vi (A) V 12 (A) 

1 

Vi 0 

Vi 2,0 

+ a( 

'Vi,/ 

W 


Vi.o 

Vj 2 ,0 

h 

.vS(A) y 2 (A). 


.^12,0 

V 2 , 0 . 

. Vjl,/ 

Vi,/ . 

. 

nlo 

Vi.o . 

\) 


(F.7) 


(F.8) 


(F.9) 
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Below, we use the notation 


Note that 


, v * dM 

M = ~dx- 


(F.10) 

(F.ll) 

(F.12) 

(F.13) 


The derivations of the expression for and are primarily based on the application 

of the following derivative formulas. It is assumed that X is an n x m matrix. 

Derivative Formulas 



A B 1 


Af — Aq Bj 


c d\ 


Cj-Co Df 

Rl 

A m ] 


*i,/ — *i,o 

*Jj *2 


.*12,/ ” *12,0 

*• 

to 


1 

o 

[VS 


L ^12,/ ” ^12,0 


-£-XA = [A(j,:)]ro*-i 

(F.14) 

-/-AX = [,4(:,t)]coi-j 

CLjC ij 

(F.15) 

/^X t A = [A(i, :)U.J 

(F.16) 

/-AX' 1 = [A(:,j)U- t 

(F.17) 

/—AXB — A(:,i)B(j,:) 

UX{j 

(F.18) 

/-AX r B = A(:,j)B(i,:) 

(F.19) 


dA 

dOj 


dA 

dl c ,kt 

dA 

dc Cy kt 


0 0 
[C(£, :)] r0 w-fc ^ - [DCc(t, :)]row-t 

0 -[B(:,A:)] co |_/ 


(F.20) 


(F.21) 
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where 


dA 


dd, 


c,kl 


- B(:,k)C(l ,:) 0 

0 0 


(F.22) 


an( j are given respectively by (D.2.36) and (D.2.37) of Appendix D. 




dR 

d9j 


dR 

db c ,kt 

dR 


= 0 


di c, 

a# 


0 [-■«»(:,*) + *)]«,!-< 


c -“ [SYM [ C c T i?2 (:^)1coi-/ + [i?2Cc (fc,:)lrow-r I 

T/. #\(d_ n r'l l . \ r?T 


dd t 


c,kl 


C D c C ( k,:)] r* ( L 

■f[CTpTn,(:,fc)-i| la (:,fc)|C«,:) ° U*J*2k£l fc > -J 


SYM 


0 


dV 

dOj 


(F.23) 

(F.24) 

(F.25) 


dV 


db 


0 


\V U {-.,1) - BDc 


c '“ [SYM [V a J? (/,:)1row-fc + \B c V 2 (:,l)Ui-k I 

av 


dc c ,k( 

dV 


= 0 


r B(:,fc)[V 3 pJfl T (<,:)-V 1 i(r,:)) 


dd, 


c,kl 


SYM 0 


j a a a 

a 1M 


A = 


A - BD c C - BD c C -BC c 
B c C -B c DC c 


i*?? 


m. 


it! = 


Rn R\2 

lT lT 

|_i? 12 ii 2 2 


(F.26) 

(F.27) 

(F.28) 


(F.29) 


(F.30) 
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where 


where 


Rn =Ri- C T DjRu ~ C T DjRj 2 - R u D c C - R u D c C 


+ C r D?R 2 D c C + C T DjR 2 D c C + C T DjR 2 D c C 

(F.31a) 

R\ 2 — — R\ 2 C c + R 2 C c ■+■ D^R 2 C c 

(F.316) 

R 22 — Cj R 2 C c - 

(F.31c) 



1 - ax 


ir = \*n V 
v 12 v 22 . 

(F.32) 

Vn = Vi- BDcVS -BD c v£ - V n Dj B t - V u DjB T 


+ BD C V 2 D C B T + BD c V 2 D c B t + BD c V 2 D c B t 

(F.33a) 

Pn = -V u Bj - BDjV 2 Bj - BD c V 2 Bj 

(F.33i») 

V 22 = b c v 2 bJ. 

(F.33c) 
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Appendix G: Derivation of and yj* * for the Input-Normal Riccati Basis 


G.l Problem Statement 

In this Appendix we assume that the controller triple ( A cy B cy Cc ) is in the input normalized 
Riccati basis described in Section 3, such that 

A e = \[C?C c - B c Bj] * I nc + [(CjC c ) * (N nc + ft T ) - (N nc + ft) * ( B c Bj ) * ft T ] *H (G.1.1) 


where 


and 


ft = diag(C* C c ) diag(£ c J5p ) 


T\-l 


H = (N nc -I n J/(ft T -ft + 7 n J. 

We derive explicit expressions for the derivatives ^ b A ‘ t and g e A c i t ■ 

Below, we use the notation 

Fb c = B c Bj 
F Cc = CjC c 
M = N nc - 7„ e 

= [1,1, — , 1] T , U nc € ffi" c 

and recognize that 

Z = diag(Z)u^ c , Z € IR ncXne 

N nc = Un c «£ 

r(k) _ e (fc) u T 

^n c e n c u n c 

ft = G/F 

H = M/(ft T - ft + 7„ c ) 

A c = ^7 * [F Ce - FbJ + [F Ce * ( N nc + ft T ) - (N nc + ft) * F Bc * ft T ] * H 
[Fc c * ( N nc + ft T ) - (N nc + ft) * F Bc * ft T ] = A e * M/H. 

The derivations of the expressions for Q b A * t and Q A t t t use the following identities. 

F (k,t) _ Jk)M) T 
J ^m x n “ c m c n 


(G.1.2) 

(G.1.3) 


(G.1.4) 
(G.1.5) 
(G.1.6) 
(G.l. 7) 


(G.1.8) 
(G.1.9a) 
(G.1.96) 
(G.l. 10) 
(G.l. 11) 

(G.l. 12) 
(G.l. 13) 
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e£ )T Z = Z(k,:), Z £ IR mXn 

(G.1.15) 

Zetf = Z(:,£), Z € IR mXn 

(G.1.16) 

±[N/D(x)] = -[N/D(x)/D(x)} * ^ 

(G.1.17) 


(G.1.18) 

Z = ab t ^ Z = (a*6)u n , a, 6, € IR-" 

(G.1.19) 

M * M = M 

(G.1.20) 
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G.2 Expression for and 

Differentiating (G.1.12) with respect to b C]k t gives respectively 
dA c 1 r dF 


dA c 1 , dF 

db cM ~ 2 n ‘ * db cM 

t „ dn T an 

+ Kc e * o7 

ubc'ki ob C 'kt 

3ff 

+ (A c * M/H) * — 

ob c ,kt 


*F Bc *n T -(N nc + n)*(^*n T + F Be *^)]*H 

Ob C 'kt ' ' 


an T , 

db c ,kt J 


(G.2.1) 


BAc _ 1 , dF^_ 
dc — 2 * ^ ^ ^ 

+ * (JVn. + fl T) + Fa . . F „. . J 5 T - (jv n , + 13 ) . ft. . j^-] . ff 

OC c ,Jt/ OC Ct kl OCc,kt OC Ci kl 


+ (4 C * M/H) * 


dc c .kl 


(GJ2.2) 


Below, we develop explicit expressions for the derivatives on the right hand sides of (G.2.1) and 
(G.2.2). 

dF Bc , dF Cc 

— and 

dbc^kt QCc^kt 

Differentiating (G.1.4) and (G.1.5) respectively with respect to b c<k e and c k< t gives 


dF Be dT i d 

K7,~ ‘ + ‘ 

r 1 i /^T 

q ~ L n u Xn c b c T 

uC Ct k£ 

which using (G.1.14) are equivalent to 

- e (*) e (') T fi T + B e ( V* )T 

Mc.fc/ n ‘ n » c c "■ n ‘ 

_ e ('U*) T <7 + C T e (fc) e ( '> T 

a. , “ C n. e n. W + W e n„ e n c • 

From (G.2.5) and (G.2.6) with (G.15) and (G.16) we obtain 


(G.2.3) 

(G.2.4) 

(G.2.5) 

(G.2.6) 


d I B ‘ = e<«ft(:,/)' r + lU-.Mf 
vvcM 

(G.2.7) 

= eWCe(fc,:) + C c (fc,:) T 4? T 

(G.2.8) 
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dQ. , dQ 

dbcM ^ fccM 

Differentiating (C.1.10) with respect to b c<k t and c cM using respectively (C.1.17) and (C.1.18) 
gives 


(Fc.IFb.IFb.)' f b B : 

UO c% kt 0”c,k t 

(G.2.9) 

an 

dc c ,kt dc c ,kt ' 

(G.2.10) 

Also, it follows form (G.2.7) and (G.2.8) with (G.1.9b) that 


d ^ Bc - 2b »e<‘)u T 
duett 

(G.2.11) 

f C * = 

OC c ,kt 

(G.2.12) 

or, equivalently, 


_ 21 £<*) 
— ZU Ct kiC'n 

OU C} k£ 

(G.2.13) 

dF Cc n c (o 

~ = 2c c,kt£n c ' 

OC c ,k£ 

(G.2.14) 

Substituting (G.2.13) into (G.2.9) gives 


^ = -(n/F Bc )*(2b cM €i l k J) 

ob c ,kt 

(G.2.15) 

or, equivalently, 


dQ _ 2 b c ,ktuk c (k) 

db c< ki Jb c ,kk nc 

(G.2.16) 


Substituting (G.2.14) into (G.2.10) gives 


an 

dc c ,kt 


2 e eM €iS/F Bc 


(G.2.17) 


or, equivalently, 





dQ 

2 C c,kt 

(G.2.18) 

dcc.kt 

1 
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dH dH 

tocM dc eM 

Let x = b C 'ki or c Ci ki- Differentiating (G.1.11) by x yields 

a -S- = -[m/(j5 t -n+ i„')heF - + /„.)]. 

It follows from (G.1.11) and (G.1.20) that 

Hence, from (G.2.20) with (G.2.16) and (G.2.18) it follows that 


dH_ m , (f W _ C<*> T ), 

OOc } kt JB c ,kk 

OCc'jd JB Cy U 


Substituting (G.2.7), (G.2.16) and (G.2.21) into (G.2.1) gives 

= -\ln. * |e<„* ) B c (:,0 T + B'(:,e)eif 

OO c k£ l 


C,Kl — 

+ [Fc. • (=7^4 ; ,T ) + ♦ -fs, . 

J B c ,kk JBcykk 

- (N nc + Q) * (e^B c (:,e) T + B^t^f ) * ti T 

- (N nc + Q) * F Bc * (-- b fb^elf)\ * H 

JBcM 

+ (A e * M/H) * [ ~ 2 ^ k t --H * H * (£<*> - 4! )T )1 

JB c ,kk 


dAc _ p (k,k) 

db c>ke cM n ‘ Xn ‘ 

+ * [F Bc * /2 T * H - Ac * M * H] 

JB cy kk 


JB cy kk 

+ 2b f’ keUk £(y * [_f c * + (JV no + J?) * F Bc * H + A c * M * H] 

JBc.kk 

rn 


- [( N nc + f2)*H*fi T ]* [e ( n k }B c (:,i) J + B c (:,ty 
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Substituting (G.2.8), (G.2.18) and (G.2.22) into (G.2.2) gives 
JET 1 - = + 

oc c ,kt 2 

+ [(<4',’C e (*, :) + C c (t, :) T eif . («„. + S 1 ) 

+ Fc. * * » 

JBc,tt JB C ,U 

- 0 N nc + Q) * F Bc * * H 

J B c ,tt 

(A, . M/H) , . H . (£<„;» - r<'/ 1 ')] 

iB c M 

or, equivalently, 

dA< _ ( e ,t ) 

rx “ C c t kt^n c xn c 
OCcM 

+ T^WA-Fbc *n T *H + A c *M*H ) 

JB Cy lt 

+ *H -(N nc + fi)*F Bc *H + A C *M*H) 

J B C ytt 

+ (JV». + &) * H * [e{‘}C c (k , :) + C c (k , :) T e£> T ]. 


(G.2.25) 


(G.2.26) 


Now, define 


H row ± [( F Bc * fi T ) - (A c * M)] * H (G.2.27) 

Hcoi = [((N nc + fl) * F Bc ) + (j4 c * M) - F Cc ] * H (G.2.28) 

H Bc = (N nc + Q)*n r *H (G.2.29) 

H Cc =(N nc + f2 T )*H. (G.2.30) 


Then, it follows from (G.2.24) and (G.2.26) that 


= -b', kl El k * l, + „(*,:) + 

OU c<k t J Bc ,kk 


dAc 

dc Ci kt 




+ tf Cc * [< } Cc(fc, :) + Cc(k, :) T e£> ] 


(G.2.31) 

(G.2.32) 
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Note that only has nonzero entries in the k th row and column, while 

entries in the / th row and column. 


only has nonzero 
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Abstract 

The minimal dimension of a linear-quadratic-gaussian (LQG) compensator is almost always 
equal to the dimension of the design plant. This deficiency can lead to implementation problems 
when considering control design for high-order systems such as flexible structures and has led to 
the development of methodologies for the design of optimal (or near optimal) controllers whose 
dimension is less than that of the design plant. This paper presents a new’ (gradient-based) ho- 
motopy algorithm for the design of reduced-order, H 2 optimal controllers. An important result is 
the development of an efficient methodology for computation of the of the cost functional Hessian 
which is required by the algorithm. The optimal controller is represented by a parameter vector 
and various parameterizations of the optimal controller are considered to reduce the algorithm 
dimensionality. The algorithm has been implemented in MATLAB and the results are illustrated 
using a benchmark, non-collocated flexible structure control problem. It is seen that the choice of a 
particular controller parameterization often introduces numerical ill-conditioning in the algorithm 
implementation. 

"Supported in part by SANDIA National Laboratories under contract 54-7C09 and NASA under 
contract NAS8-38575. 




1. Introduction and Nomenclature 


The linear-quadratic-gaussian (LQG) compensator (Atlians 1971, Kwakernaak and Sivan 1972) 
has been developed to facilitate the design of control laws for multi-input multi-output (MIMO) 
systems. An LQG compensator minimizes a quadratic performance index and (under mild condi- 
tions) is guaranteed to yield an internally stable closed-loop system. Unfortunately, however, the 
minimal dimension of an LQG compensator is almost always equal to the dimension of the plant 
and can thus often violate practical implementation constraints on controller order. This deficiency 
is especially highlighted when considering control- design for high-order systems such as flexible 
space structures. Hence, a very relevant area of research is the development of methodologies that 
will enable the design of optimal controllers whose dimension is less than that of the design plant 
(i.e., reduced-order controllers). 

Two main approaches have been developed to tackle the reduced-order design problem. The 
first approach attempts to develop approximations to the optimal reduced-order controller by reduc- 
ing the dimension of an LQG controller (Yousuff and Skelton 1984a, Yousuff and Skelton 1984b, 
Anderson and Liu 1989, De Villemagne and Skelton 1988, Liu, Anderson and Ly 1990). These 
methods are attractive because they require relatively little computation and should be used if 
possible. Unfortunately, they tend to yield controllers that either destabilize the system or have 
poor performance as the requested controller dimension is decreased and/or the requested authority 
level is increased. Hence, if used in isolation, these methods do not yield a reliable methodology 
for reduced-order design. 

The second approach attempts to directly synthesize an optimal, reduced-order controller by 
a numerical optimization scheme (Levine, Johnson and Athans 1971, Martin and Bryson 1980, 
Mukhopadhyay, Newsom and Abel 1982, Ly, Brysib and Cannon 1985, Mukhopadhyay 1987, Kuhn 
and Schmidt 1987. Richter 1987, Makila and Toivonen 1987, Kramer and Calise 1988, Mukhopad- 
hyay 1989, Richter and Collins 1989, Mercadal 1991 ). Almost all of these schemes are gradient-based 
parameter optimization approaches; that is, they represent the controller by some parameter vector 
and attempt to find a vector for which the gradient of the cost functional is zero. 

With the exception of Mercadal 1991, all of the previous, gradient-based optimization tech- 
niques are descent methods. That is, at each iteration the cost function is decreased. An alternative 
(Mercadal 1991) is to develop a gradient-based homotopy algorithm that allows an initial controller 
to be deformed gradually into a desired optimal controller by following a homotopy path. This type 
of algorithm is distinct from the previous algorithms in that each iteration does not necessarily de- 
crease the cost function. In fact, the cost may actually increase as the homotopy path is traversed. 
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However, it is quite possible that the shortest path from the initial controller to the desired optimal 
controller is not a descent path. 

Efficient path following requires accurate computation of the Hessian of the cost functional. 
Hence, this paper develops an efficient method for computing the Hessian. An alternative method 
for computing the Hessian is presented in an earlier publication (Sun 1991). However, to our 
knowledge, this previous method, based on the results of Sun 1990, does not exploit certain low 
rank matrices as does the method presented here. 

A homotopy algorithm for optimal reduced-order design is described in Richter 1987 and Richter 
and Collins 1989. This algorithm is based on solving a set of “optimal projection” equations (Hyland 
and Bernstein 1984, Haddad 1987) that are a characterization of the necessary conditions for 
optimal reduced-order control. Unfortunately, the algorithm has sublinear convergence properties 
and the convergence slows at higher control authority levels and may fail. Homotopy algorithms 
for optimal reduced-order modeling, based on optimal projection equations, are discussed in Zigic 
et al. 1992 and Zigic et al. 1993. These algorithms are based on more efficient path following 
techniques but are relatively slow due to the large dimensionality of the algorithm formulation. 

This paper describes a homotopy algorithm for the design of reduced-order, Ho optimal con- 
trollers which is not based on the optima] projection equations. The algorithm relies on the first 
and second derivatives (i.e., the gradient and Hessian) of the cost functional with respect to a pa- 
rameter vector describing the controller and an efficient methodology for computing the Hessian is 
developed. To reduce the dimensionality of the algorithm, various parameterizations of the optimal 
controller are considered. The algorithm has the potential for quadratic convergence rates along 
the homotopy curve. The results have been implemented in MATLAB and are illustrated using 
a benchmark, non- collocated flexible structure control problem. It is seen that the choice of a 
particular controller parameterization often introduces numerical ill-conditioning in the algorithm 
implementation. The algorithm presented here is similar to that described in Mercadal 1991. How- 
ever, whereas Mercadal 1991 focuses on theoretical issues related to homotopies and only describes 
a rudimentary homotopy algorithm, the present paper focusses on numerical algorithmic issues and 
describes a much more refined and efficient homotopy algorithm. 

The paper is organized as follows. Section 2 describes the Hj optimal reduced order dynamic 
compensation problem. Section 3 gives a brief overview of homotopy methods. Section 4 then 
develops a homotopy algorithm for the design of reduced-order Ho optimal controllers. Section 5 
applies the algorithm to a benchmark structural control problem and compares various algorithm 
options. Finally, Section 6 presents the conclusions. 
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Nomenclature 


Y > Z 

Y — Z is nonnegative definite 

Y > Z 

Y — Z is positive definite 

z ijiZi,j or 

(i,j) element of matrix Z 

It 

r x r identity matrix 

tr Z 

trace of square matrix Z 

vec(-) 

the invertible linear operator defined such that 

vec(s) = [s?" sj • • • s*] T , S 6 nt 7 ’* 9 
where s : £ IR P denotes the j th column of S. 

Cm 

the m-dimensional column vector whose t lh element 
equals one and whose additional elements are zeros. 

E (k<) 

-*- / m ,n 

the m x n matrix whose (k y t) element equals one 
and whose additional elements are zeros. 

z(t,o 

k ih row of the matrix Z 
(MATLAB notation) 

Z(:,i) 

k ih column of the matrix Z 
{MATLAB notation) 


2. H 2 Optimal Reduced-Order Dynamic Compensation 

Consider the system 

x(0 = Ax(t) + Bu(t) + ^(t) (2.1c) 

2/(0 = Cx(1) + Du{t) + w 2 (t) (2.16) 

where x £ IR n *,u £ IR n * , > 2/ G TR n *,wi £ IR n ' is white disturbance noise with intensity V\ > 
0, w 2 € DC 1 * is white observation noise with intensity V 2 > 0, and icj and xv 2 have cross correlation 
Vi 2 € IR n ** n, \ We desire to design a fixed-order dynamic compensator. 


i c (*) = A c z c (t) + B c y{t) 

(2.2a) 

v(t) = -C c x c (t) 

(2-2 b) 


which minimizes the steady-stale performance criterion 

J(A c ,B c ,C c )= Urn E[i T (t)R l T({) + 2T T (t)R 12U (l) + v r (t)R 2 u(t)] (2.3) 

t— OC * 
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where x c G IR nc ,n e < n,R 2 = Rj > 0 and Rj — R . J > 0. We will call this problem the optimal 
reduced-order dynamic comjtensation problem. 

The closed-loop system corresponding to (2.1) and (2.2) can be expressed as 


i(t) = Ax(t ) + w(t) 


(2.4) 


where 


x(t) = 


x c (t) 


. «H0 = 


Wy(t) 

[B c tn 2 (t) 


A* 


A 

[B c C A c 
In addition, the cost (2.3) can be expressed as 


-BC c 
- B c DC c . ' 


J(A Cl B c ,C c ) = lim £[x T (0J?i(0] 

f— *oo 


where 


R± 


R i 

CjRl 


RnC e 1 
J 2 CjR 2 C c J* 


(2.5, 6, 7) 

( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 


To guarantee that the cost J is finite and independent of initial conditions we restrict our 
attention to the set of stabilizing compensators, S c = {(A c , B c ,C c ): A is asymptotically stable}. 
Assume (A c , B c ,C c ) G S c and define Q G IR nxXn * and P G Dl" xX " x to be the closed-loop steady- 
state covariance and its dual, i.e., 

0 = AQ + QA T + V (2.11) 

0 = A t P + PA + R 


where 


Then, the cost can be expressed as 


y± \ VnBj ] 

~ [b c v? 2 b c v 2 bJ j- 


J{A c ,B c ,C c )= tr QR = tr PV. 


Also, Q and P can be expressed in the partitioned forms 


( 2 . 12 ) 

(2.13) 

(2.14) 


Q = 


5u Qn 

Vn O22 J 


<?„ g nt nxXnx , (? 22 e m. n ‘ xn ‘ 


p = 


[JV p 2 \ > Aiffi n * xn ‘, 


(2.15) 

(2.16) 
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Notice that Qn is the covariance of the plant states, Q 22 is the covariance of the compensator 
states and Qu is the cross- covariance of the plant and controller states. 


Expressions for the partial derivatives and are given below. First, we define Z 

satisfying 

Z = QP (2.17) 

and assign Z the partitioned form 


Z = 



Z\2 

Z22 ' 


Z„eIR n ' xn *, Z 22 € IR n ‘ Xn ‘, 


(2.18) 


The cost derivatives are then given by 


dJ 

dA e 



(2.19) 


jg- e = 2(Fj 2 Vj2 + P 22 B e V 2 + Z 1 T ,C T - Zj 7 CjD T ) (2.20) 

^ = 2(-Rj 7 Q l7 + RiC c Q 22 - B t Z 2 t , - D T BJZJ 2 ) (2.21) 


Definition 2.1. A compensator (A c , B c ,C c ) is an extremal of the optimal reduced-order 
dynamic compensation problem if it satisfies the stationary conditions 





( 2 . 22 ) 


The homotopy algorithm of Section 4 is based on finding extremals of the optima] reduced-order 
dynamic compensation problem. 


3. Homotopy Methods for the Solution of Nonlinear Algebraic Equations 

A “homotopy” is a continuous deformation of one function into another. Over the past several 
years, homotopy or continuation methods (whose mathematical basis is algebraic topology and 
differential topology (Lloyd 1978)) have received significant attention in the mathematics litera- 
ture and have been applied successfully to several important problems (Avila 1874, Wacker 1978, 
Alexander and Yorke 1978, Garcia and Zangwill 1981, Eaves, Gould, Peoitgen, and Todd 1983, 
Watson, 198G). Recently, the engineering literature has also begun to recognize the utility of these 
methods for engineering applications (see e.g. Richter and DeCarlo 1983, Richter and DeCarlo 

1984, Turner and Chun 1984, Dunyak, Junkins, and Watson 1984, Lefebvre, Richter and DeCarlo 

1985, Sebok, Richter, and Decarlo 198G, Horta, Juang and Junkins 1986, Kabamba, Longman and 
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Jian-Guo 1987, Shin, Raftka, Watson, and Plautt 1988, Rakowska, Haftka, and Watson 1991)). 
The purpose of this section is to provide a very brief description of homotopy methods for finding 
the solutions of nonlinear algebraic equations. The reader is referred to (Watson 1986, Richter and 


DeCarlo 1983, Watson 1987, Watson 1986) for additional details. 

The basic problem is as follows. Given sets 0 and $ contained in IR n and a mapping F : 0 — * $, 
find solutions to 

F{6) = 0. (3.1) 

Homotopy methods embed the problem (3.1) in a larger problem. In particular let H : 0 x [0, 1] — ► 
IR n be such that: 

1) ff(0,l) = F(6). (3.2) 

2) There exists at least one known 6 0 G Dt" which is a solution to H(,0) = 0, i.e., 

H(0 o ,0) = 0. (3.3) 

3) There exists a continuous curve (0(A), A) in IR" x (0,1) such that 

/f (0(A), A) = 0 for A € [0, 1] (3.4) 

with 

(*(O),O) = (0o,O). (3.5) 


4) The space 0 X [0, 1] has a differential structure so that the curve (0(A), A) is differentiable. 

A homotopy algorithm then constructs a procedure to compute the actual curve a such that the 
initial solution 0(0) is transformed to a desired solution 0(1) satisfying 

0= tf(0(l),l) = -F(0(1)). (3.6) 


Differentiating H(6( A), A) = 0 with respect to A yields Davidenko's differential equation 


OH d_B OH 
06 dX + OX 


(3.7) 


Together with 0(0) = 0o, (3.7) defines an initial value problem which by numerical integration from 
0 to 1 yields the desired solution 0( 1 ). Some numerical integration schemes are described in Watson 
1986 and Watson 1987). 


C 



4. A Homotopy Algorithm for H 2 Optimal Reduced-Order Dynamic 
Compensation 


This section presents a new homotopy algorithm that can be used to design H 2 optimal reduced- 
order dynamic compensators. Particular attention is given to construction of the Jacobian of the 
homotopy map. 


9* 


4.1 The Homotopy Map 
If we define 

vec(4 c )' 

vec(5 c ) I , (4.1) 

vec(C c ) J 

then the cost functional of Section 2 can be expressed as J(9). The homotopy defined in this section 
is based on finding 9 satisfying 

. ft J 

(4.2) 


0 = m ^ %m- 


It is useful to recognize that 


d£ 

06 

dj dJ 


vec 


vec 


vec 


QJ n 
dA* 

JLL 

dB c 

8J 

dC c J 
8J 


(4.3) 


Expressions for the partial derivatives and are given by (2.19)-(? 21). 


Definition of the homotopy map H(6 y A) 


To define the homotopy map we assume that the plant matrices (4, B, C, D ), the cost weighting 
matrices (J?i, #2* # 12 ) and the disturbance matrices {V\yV?,V\ 2 ) are functions of the homotopy 
parameter A G [0, 1]. In particular, it is assumed that 


A(A) 5(A) 
C( A) D( A)J 


Ao 

Co 




Ao Bq \ 

Co D 0 \J ’ 


(4.4) 



5i(A) 5i 2 (A) _ f /nrT/n 

R? 2 (\) 5o(A)J - Lh(A)L « (A) 

(4.5) 

where 

Lr( A) = Lftfl 

+ ^{Lrj - Lfifi) 

(4.6) 

and Lfi t o and Lrj satisfy 


5 10 #12,0 ] 

. ^12 ,0 ^ 2.0 J 

(4.7) 


II 

#1,/ #12, /I 

.^12./ ^2./ J 

(4.8) 
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v'{(\) V?\\) - Lv(*)Ll(\) 


(4.9) 


where 





Lv( A) = Lv,o + X{Ly,j - Lv, 0 ) 

(4.10) 

and Lv, 0 and Ly,j satisfy 

Lv,oLy 0 — 

■ V li0 1*12,0 1 
.vg, 0 V 2 . 0 J 

(4.11) 


LvjLyj = 

Vi.f Vn,f] 
[Vil / ^,/J- 

(4.12) 


Note that (4.4)-(4.12) imply that /4(0) = Aq and >1(1) = Aj , 2?(0) = Bq and f?(l) = Bf, etc ... 

and it is understood that Aj,Bj ,. . . were referred to previously simply as A,B, The change in 

notation is simply for convenience. 

Let P( A), <3(A) and Z{ A) satisfy 

0 = A(\) T P(\) + P(A)i(A) + R( A) (4.13) 

0 = i(A)Q(A) + <?(A)y4(A) T + V(A) (4.14) 

2(A) = <?(A)P(A) (4.15) 


with partitioned forms 


P( A) = 


PuW 

A t 2 (a) 


Pn( A)' 
p 22 (a)J’ 


Q( A) = 


Q n(A) 


Ql2(A)' 
<? 22 (A) ’ 



2n(A) 

222 (A) 


(4.16) 


where the (1,1) and (2,2) blocks of each matrix are respectively n x x n x and n c x n c . The homotopy 
map H(8, A) is defined as the gradient of the cost of the system defined by the homotopy parameter 


A. In particular, 


where 


B {8, A) = 


-vec^OU))' 

vec(P Bc (#,A)) 

.vec(P Ce (0,A)). 


(4.17) 


B At (S,X) = 2Zj 2 

BbA^X) = 2(P, t 2 V 12 + P 22 B C V 2 + Zj 2 C T - Zj 2 CjD T ) 
Bc c (0, A) = 2(-P? 2 <? 1 2 + P 2 C c Q2 2 - B T Z 2l - D T BjZ? 2 ) 


Note that in (4.18)-(4.20) and below the argument A is omitted for notational convenience. 


(4.18) 

(4.19) 

(4.20) 



4.2 The Jacobian of the Homotopy Map 


We now consider that computation of the Jacobian of H(9, A). Note that 



VH(6,\) r = [He Hx) 

(4.21) 

where 

„ * dH u * OH 

Ht = Je' Hx = J\- 

(4.22) 

Since H(0 , A) is the gradient for the system defined by A, H$ is the corresponding Hessian. Recalling 

that 6 is defined by (4.1), such that for some integers k and i, 8, is given by 



Qj ~~ a c ,ki> 0j — ,kh &j — C c t kh or 0j “* 

(4.23) 

It follows from (4.13) that H$ is of the form 


He = 

■ • • vec ( a.?.. , Ha. ) • • • vec( S6 f 4 , H Ac ) . . . vec( 8c f w H ) 

• • ' vec ( H B < ) • • • vec ( e*J M vec ( 9c a c ki H Bc ) 
L- vec ( *<?,)•.. v ec(a^/f Ct ) • • • vec( 8c f w H Cc )\ 

(4.24) 

and H \ can be expressed as 

' ve c(jx H *<Y 

H x = vec (£ x H B c ) . 

L vec ( "ax ) . 

Below, we develop explicit expression for the derivative terms appearing on the 

(4.25) 

right hand sides of 

(4.20) and (4.21). We use the notation 



M U) A 

09, 

M * 8M 

Um W 

(4.26) 

(4.27) 

Differentiating (4.13)— (4.15) with respect to 0 : yields 



0 = A T P U) + P U) A + (A (j)T P + PA U) + R U) ) 

(4.28) 


0 = AQ (}) + Q^A 1 + (A U) Q + QA (i)T + V U) ) 

(4.29) 

Z U) = Q^P + QP U) 

(4.30) 


where expressions for the derivatives A^\R.^ and are given by (A.20)-(A.28) of Appendix 
A. Similarly, differentiating (4 . 13 )— (4. 15) with respect to A yields 

o = a t P + Pa + (a T P + PA + P) 

0 = A$ + ClA 1 + (AQ + QA + &) 

Z - QP + QP 
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(4.31) 

(4.32) 

(4.33) 



where expressions for A y R and V are given by (A.29)-(A.33) of Appendix A. 

Before presenting the desired derivative expressions we define 

H‘ At (Zl»)±2Zif (4.34) 

H' Bc (P u \Z u) ) = 2(Pjf T Via + P^BcVt + C* 1 - tyfcjD*) (4.35) 

ff'c c (Q U) J U) ) = 2 (-Rj t Q l £ + RiCcQg - B T Z<f - D T BjZ$f) (4.36) 


Notice that the right hand sides of (4.34)-(4.36) are identical in form to the right hand sides 
of (4.18)-(4.20). The only difference is that P,Q, and Z have been replaced by and Z M 

Derivatives with respect to a c \ki 


Differentiating (4.18)-(4.20) with respect to a c ,fcf( = ffj) gives 


OHa « 

da.c,kt 

9Hbc_ 

da c ,kt 

dHc, 

da Cikt 


H'aAZ U) ) 

H' Bt {P u \Z (i) ) 


(4.37) 

(4.38) 

(4.39) 


Derivatives with respect to b c ,kt 
Differentiating (4.18)-(4.20) with respect to b c ,ki(— 0j) gives 


- H ^ zh)) 

(4.40) 

= H' Bc (P^\Z^) + 2P n E[ k /J n V 2 

(4.41) 

= H'cSQ (i \Z™) - 2 D T EH; k J nc zJ 2 . 

(4.42) 


Derivatives with Respect to c Ci kt 

Differentiating (4.18)-(4.20) with respect to c c ,*/(= 8j) gives the following. 

= H\ c {Z (3) ) 


9JU rr, 

dc c ,kt 
Mb. 

dc c . 


kt 


= H' B (P U) ,Z U) ) - 2Z, T 2 4!xl. D 1 


= H'cSQ u) ,ZM) + 2B 2 Eiy nt Q 22 . 


(4.43) 

(4.44) 

(4.45) 
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Derivatives with respect to A 


Differentiating (4.18)-(4.20) with respect to A gives 


dH Ac 

dX 

dH Bc 

dx 

dHcc 

dX 


= H'aSZ) 

= H' Bc (Pj) + 2(P?V U + P 22 b c v 2 + z£c T - Zj 2 Cji ) r ) 

= H' Cc (ti, ’z) + 2(-PJ 2 Q u + R 2 C c Q 22 + B T Z7 2 - i> r BjZj 2 ) 


where from (4.4)-(4.12) 


A B 
C D 


R\ R\i 

M 


A j — Aq Bj — Bq 
C/ — Co Cf — Co 

= LrLf? + LrLr 


where 


where 


Lr — Lrj — IRQ 


V } V u 

'ft v. 


= Ly Ly + LyL\ 


Ly = Lyj — Ly % o« 


( 4 . 46 ) 

( 4 . 47 ) 

( 4 . 48 ) 

( 4 . 49 ) 

( 4 . 50 ) 

( 4 . 51 ) 

( 4 . 52 ) 

( 4 . 53 ) 


He can now be computed using (4.24) and (4.37)-(4.45). 

Note that the calculation of the j th column of He requires the computation of the Lyapunov 
equations described by (4.28) and (4.29). Significant computational savings can be made by solving 
these Lyapunov equations in a basis in which the closed-loop state matrix A is nearly diagonal (i.e., 
a modal form) or nearly block triangular (i.e., a Schur form). This requires transforming the 
corresponding forcing terms into this basis which can be costly if the dimension of the closed-loop 
system, n c /(= n x + n c ) is large. In fact, if the forcing terms are dense, this transformation requires 
2 operations. Fortunately, it is seen by (A.20)-(A.28) of Appendix A that these forcing terms 
are low rank. Hence, these transformations do not have to be expensive and often require only 
about 2 nit operations. Computation of the expressions (4.37)-(4.45) requires the solutions of the 
Lyapunov equations in their original basis. However, it is not efficient to numerically perform this 
transformation before substituting into (4.37)-(4.45). Instead, symbolic substitution and judicious 
choice of the order of matrix multiplications can result in significant computational savings. The 
details of efficient computation of He are presented in Appendix B. 
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H\ is computed using (4.25) and (4.46)-(4.48). This requires computation of the Lyapunov 
equations (4.31) and (4.32). The forcing terms for these Lyapunov equations are not sparse so 
that computing H\ in a particular basis requires 2n^ operations to transform the forcing terms. 
However, the rest of the optimization associated with the computation of H$ does apply to the 
computation of H \ . 

4.3 Reduction of the Dimension of the Controller Parameter Vector (9) 

The homotopy function H (0, A), described earlier, was defined to solve the H 2 optimal reduced- 
order dynamic compensation problem. The vector 9 was defined such that it contained each of the 
elements of the controller matrices, A c , B c and C c . However, for computational efficiency it is 
desired that 9 be as small as possible. Hence, we desire to represent the controller matrix with the 
fewest parameters possible (i.e., we desire 9 to have the smallest dimension possible). The minimal 
number of parameters p, n i n with which a compensator can be represented is given by (Martin and 
Bryson 1980, Denery 1971) 

Pmin = n c (n u + n v ) (4.54) 


One canonical form which allows representation of a controller with a minimal number of 
parameters is the modal form described in (Martin and Bryson 1980). This form will be called 
here the Second-Order Polynomial (SP) form. For this parameterization a triple (A c . B c ,C c ) has 
the following structure. 

A c = block- diag{y4 Cili y4 c . 2 ...,/l Cir } (4.55) 


where A C: i £ IR 2x2 for i € {1,2, ...,r} and each A c ,i (with the exception of A CiT if the row 
dimension of A c is odd) has the form 


A 


c,i 


o 

o> 


L "c.i 



(4.56) 


to allow for either a complex conjugate set of poles or two real poles. B c is completely full and 


Cc — [Cc,„ C Ci 2 , • - . , C c ,r] 


where C c ,i has the form 


Cc. r = 


1 0 
* * 

♦ * 


(4.60) 


(4.57) 
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The controller canonical form described in Kailath 1980 also allows representation of a controller 
with a minimal number of parameters. For single-input, single-output (SISO) systems in controller 
canonical form the A c matrix is a companion matrix. In particular, A c has the form 


A c = 


'0100 
0 0 10 
0 0 0 1 


O' 

0 

0 


* * * * • • • * 


(4.58) 


In. addition, 


and C c is completely full. A dual form 
form (Kailath 1980). 


'O' 

0 


B c = 


0 


(4.59) 


LiJ 

of the controller canonical form is the observable canonical 


It is also possible to represent the controller in a basis where the number of free parameters p 
satisfies 

Pmin < p < Pmax = ^c(«c + Tl u + 7l y ). (4.60) 

One such basis is the tridiagonal basis (Geist 1991, Parlett 1992) in which the controller state 
matrix is constrained to have nonzero elements only on the diagonal, the super-diagonal and the 
sub-diagonal. That is, 

r* * 

♦ * 0 

A c = | * * (4.61) 

* * - 

B c and C c are completely full. For this form the number of free parameters is given by 

V = 7*ruin “h (dn c — 2) 


A common feature of each of the above bases is that they are described by simply constraining 
certain elements of the controller (or plant) matrices to constant values (e.g., 1 or 0) while allowing 
the remaining parameters to have arbitrary values ( A c , B c , C c ). Hence, the corresponding parameter 
vector (0 p ), gradient vector and Hessian matrix (/f* )P ) are given by 


0 P = TO 

(4.62) 

Je,p = 

(4.63) 

h,, p = r#,r T 

(4.64) 
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(4.63) 

(4.64) 


Jq,p — TJ$ 

H e , p = r#,r T 

where T is an elemental matrix (i.e., each row has only one nonzero element and this element has 
unity value). It should be noted here that H$ %p can be computed more efficiently than shown in 
(4.64). Since it is not necessary to construct the large Hessian H$ to compute the smaller Hessian 

H»,p- 


4.4 Overview of the Homotopy Algorithm 

This section describes the general logic and features of the homotopy algorithm for Hi optimal 
reduced-order control. It is assumed that the designer has supplied a set of system matrices, 
Sj = ( A j , B / , C / > D j > R\ j , Ri j , Rii j j y Vi j ,Vu t j) describing the optimization problem whose 
solution is desired. In addition, it is assumed that the designer has chosen an initial set of related 
system matrices S 0 = (Ao, 01 ^ 2 , 01 ^ 12 , o^.o^.o^I^o) that has an easily obtained 

optimal controller (4 Ct o, B Ct o,C c ,o) of order n c . 

It is always possible to choose the initial system So such that (4o, Bo, Co, Do) in nonminimial 
with minimal dimension n c . In this case, it is easy to show that the corresponding LQG compensator 
has minimal dimension n r < n c and will usually have minimal dimension n r = n c . In the latter 
case, (A c ,o, B c , 0 )C c ,o> D c ,o) is chosen as a minimal realization of the LQG compensator. However, 
we have seen experimentally that the corresponding homotopy can lead to failure of the homotopy 
algorithm. Similar observations have been made by Mercadal (Mercadal 1991). In particular, 
Mercadal has shown that allowing the plant parameters to vary along the homotopy path can lead 
to the development of destabilizing controllers or path bifurcations. 

That the above type of homotopy would cause problems is somewhat intuitive since for a given 
A, say Ai 6 (0,1], a controller (A^X^). B c (Xi),C c (\\)) that stabilizes the plant (4(A X ), B(Aj ), 
C(Ai),Z?(Aj)) may not stabilize the plant (4(A 2 ), i?(A 2 ),C(A 2 ), Z)(A 2 )) for A 2 ^ A x . Hence, below 
we present ways of constructing the initial system So that does not require the plant paramaters 
(A,B,C,D) to vary along the homotopy path. In this case, a controller that stabilizes the plant at 
Ai will also stabilize the plant at A 2 > A x . This argument in itself does not ensure that at every 
step along the homotopy algorithm the controller design remains stabilizing. This is a subject that 
requires further research. It should be mentioned that another advantage of a homotopy that varies 
only the performance weights (/? x ,/? 2 , J? x2 , Vj, V 2 , V] 2 ) is that the optimal controller at each point 
is optimal with respect to the real nominal plant (4/, JB/, C/, Dj). 


14 


Now, we present three options for constructing So and hence defining the horaotopy. 

Option 1. One alternative for constructing So is to choose Ao to be stable (e.g., if A/ is stable, 
let Aq = Aj or if Aj is unstable, let Ao = Aj — a I where a is sufficiently large to ensure stability of 
Ao), and let either o or Vi f o be zero with all other parameters equal to their final values. In this 
case (A Ci o, 5 C| o,C C) o) is chosen such that it’s input-output map is zero, i.e., C Ct o(sI nc — A Ci o)“ 1 iJ*,D 
= 0 . 


Option 2. Another alternative is to choose Ao to be stable and as elaborated in (Collins, 
Haddad, and Ying 1993) and choose either ^ 2 , 0 ) or (Vi t o,i? 2 ,o) as given below. (Again, all 
other initial parameters are equal to their final values.) 


(i) In a basis in which 

Ao = 

choose R l0 to be of the form 


(Ao)n 0 
L(Ao)21 ( Aq )22 J 


1 (Ao)n € m 


n c x n e 


Rifi = 


(Ri,o)u 0 
0 0 


, (Ri,o)n e ffi. 


n c X n c 


( 4 . 65 ) 


( 4 . 66 ) 


and for some positive scalar a choose 


^2,0 = 0V2./ 


( 4 . 67 ) 


(ii) In a basis in which 


An = 


(Ao)n (Ao)i 2 

0 (Aq)22 


, (Ao)n eiR n ' XT \ 


choose V\ t o to be of the form 

V 1J0 = 

and for some positive scalar a choose 


,o )i 1 0 

0 0 


, (Vi f o)ii € IR 


n c X n c 


( 4 . 68 ) 


( 4 . 69 ) 


ft 2 — 


( 4 . 70 ) 


As discussed in (Collins, Haddad, and Ying 1993), a in (4.67) or (4.70) can always be chosen 
sufficiently large that the corresponding LQG compensator is nearly nonminimal. In this case, a 
very close approximation to (A Ct o, 5 Ci o, C Ci o) is easily obtained by reducing the LQG compensator 
to it’s (nearly) minimal realization using an appropriate technique such as balanced controller 
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reduction (Yousuff and Skelton 1984). This initialization option can sometimes present a shorter 
path to optimal solution than the first option given above. 

Option 3. A third alternative (which does not require Aq to be stable) is based on the following 
experimental observation. The initial system can be chosen to correspond to a low authority control 
problem, e.g., one can choose 

# 2,0 = ^ 2,0 = PVi j ( 4 * 71 ) 

with a and large and let all other initial system parameters equal their final values. In this case 
it has been observed that the reduced-order controller (A c>r , J9 Cf r,C Ct r) obtained by suboptimal 
reduction of an LQG controller will often yield virtually the same cost as the LQG controller (see, 
e.g., De Villemagne and Skelton 1988), hence indicating that (A Ctr ,B Cir ,C Cir ) is nearly optimal. In 
this case we choose (A Ct o,B Ct o,C Ct o) = (A c , r , B c ,r>C c ,r)- It should be noted that these observations 
are partially (but not fully) explained by the results of (Collins, Haddad, and Ying 1993). 

Below, we present an outline of the homotopy algorithm. This algorithm describes a predic- 
tor/corrector numerical integration scheme. There are several options to be chosen initially. These 
options are enumerated before presenting the actual algorithm. Notice that each option corresponds 
to a particular flag being assigned some integer value. 
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Controller Basis Options : 


basis = 0. No basis (i.e., all elements of the controller matrices are considered free.) 

basis = 1. Tridiagonal Basis. 

basis = 2. Second-Order Polynomial Form. 

basis = 3. Controllable Canonical Form. 


Note that for basis = 0 or 1, p > p m -,„ while for basis = 2 or 3, p = p, nin . 

Prediction Scheme Potions : 

Here we use the notation that Ao,A_i, and Aj represent the values of A at respectively the 
current point on the homotopy curve, the previous point and the next point. Also, 8 p = d8 p /dX 
and is the solution of Davidenko’s differential equation (4.7), rewritten here as 


H*, P 8' P W + H X = 0 . 


(4.72) 


If P = Pmim He, P is generally invertible, then 6 p ( A) is given exactly by 



(4.73) 


If P > pmin, then He <p generally has rank p m and 6 p ( A) is approximated by the least squares 
solution of (4.73) or 


e' p = -V 


Zo l 0 
0 0 


U T II X 


where it is assumed the IIe, P has the singular value decomposition 


H e , P = U 


—•o u 

0 0 


V'T, S 0 £ xp mill 


(4.74) 


(4.75) 


Note that for p = p mm (4 .73) and (4.74) are equivalent. 


pred = 0. No prediction. This option assumes that 8 p (Xi) = 0 P (Ao). 

pred = 1. Linear prediction. This option assumes predicts 0 p (Ai) using only 0 p (A o ) and 8 p '(Xo). 
In particular, 

^p(^i) = 0p(Ao) + (A) — AoJtfp^Ao) (4.76) 


pred = 2. Cubic spline prediction. This option predicts 0 p (Aj) using 8 p ( Ao), ^ p '(A 0 ), 0 p (A_j) 
and V(A-i). In particular, 

^p(Aj) = ao + «i Ai + ojA) 2 + asAj 3 (4.77) 
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where a 0 , aj, a 2 and 03 are computed by solving 


[o 0 O] a 2 a 3 ] 


■ 1 

0 

1 

0 ■ 


r«p(A-.)-| 

A_i 

i 

^0 

1 



A 2 -, 

2A_j 

^0 

2Ao 


1 

n.-. a. 
1 

Lai, 

3Al, 

Ag 

3Aq J 



(4.78) 


Note that if this option is chosen, then at the initial algorithm prediction step 0 p (A_i) and 
0 p(A-i) are not available, in which case linear prediction is used. 

Correction Options : 


Here we assume that the homotopy parameter has a fixed value Ao. The vector 0 P represents 
the current approximation of the parameter vector at A = Aq. Each of the options corresponds to 
updating 6 p using the formula 



0 p *— 0 p + Aflp 

(4.79) 

where 




A0 P = —Ge, p Je, P 

(4.80) 

for some choice of G$ tP . 




corr = 1 . Newton correction. In this option, if p = 

G$, v = (4.81) 

while if p > p min , 

Ge, P = V(E 2 + a 2 iy l ZU T (4.82) 

where a is some (small) scalar and (U, V, E) denote the singular value decomposition of He, p 
such that 

H e%P = UZV r . (4.83) 


It can be shown that if Ge >p is given by (4.82), then minimizes 

+ t r f + « ! l|Afl,|| J ). (4.84) 

Hence, is essentially a “Newton correction” that is relatively insensitive to singularities or 
near singularities in the Hessian, He, P - 

corr = 2 . Quasi-Newton correction. In this option, Ge tP denotes the estimate of using 
only gradient and cost information. For the algorithm presented here the BFGS inverse Hessian 
update is used (Fletcher 1987). 
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Outline of the Homotopy Algorithm 

Step 1. If basis > 1, then transform the initial controller (A Ct o, B c $, C C) o) to the chosen basis 
and let 6 0 ,p be the corresponding vector of free parameters. 

Step 2. Initialize loop = 0, A = 0, AA G (0,1], 5 = So, & P = 0o,p and compute the cost J 
and the cost gradient J 6 p corresponding to S and the controller described by 0 p . 

Step 3. Let loop = loop+1. If loop = 1, then go to Step 5. Else, continue. 

Step 4. Advance the homotopy parameter and predict the corresponding parameter vector 0 
as follows. 


4a. Let Ao = A 
4b. Let A = Ao *4* AA. 

4c. If pred > 1, then compute 0J,(Ao). 

4d. Predict 0 P (A) by using the option defined by pred. 

4e. If the normabzed gradient ^ lP ||C*0 lP ||/||0 p || satisfies some preassigned tolerance, then 
continue. Else, reduce AA and go to Step 4b. 

Step 5. Correct the current approximation 9 V to the optimization problem defined by S using 
the option defined by corr until the normalized gradient, 


\\M\ 

satisfies some preassigned tolerance. 

Step 6. If A = 1, then stop. Else, go to Step 3. 


(4.85) 


The above algorithm assumes monotonicity of the solution curve as a function of the homotopy 
parameter A. However, it is not difficult to modify the algorithm so that the variable parameter is 
the arc length as discussed in Watson 198G and Watson 1987 since this modification would still only 
require the computation of H$ and H\. The modified algorithm would not require monotonicity of 
the solution curve. However, so far in our computational experience the solution curve has always 
been monotonic. 

Note that if p = p min and corr = 1, then the corrections of Step 5 correspond to Newton 
corrections. Hence if the prediction tolerance used in Step 4 is sufficiently small, then, entering 
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Step 5, 6 p will be close enough to the optimal value 6* so that the quadratic convergence proper- 
ties of Newton’s method (Fletcher 1987) can be realized. In practice, this quadratic convergence 
property is not always realized due to numerical ill-conditioning associated with the minimal pa- 
rameterization of the controller. This ill-conditioning is illustrated and discussed further below. 

5. Illustration of Reduced-Order Design Using a Four Disk Example 

This section illustrates the homotopy algorithm of Section 5 by considering control design for 
an axial beam with four disks attached as shown in Figure 5.L This example was derived from 
a laboratory experiment described in (Cannon and Rosenthal 1984) and has been considered in 
several subsequent publications [Anderson and Liu 1989, De Villemagne and Skelton 1988, Liu, 
Anderson, and Ly 1990, Hyland and Richter 1990). The basic control objective for the four-disk 
problem is to control the angular displacement at the location of disk 1 using a torque input at the 
location of disk 3. It is also assumed that a torque disturbance enters the system at the location 
of disk 3. 

The design philosophy adopted here is that the scaling q2 of the nominal control weight i? 2 ,o = 1 
and the nominal sensor noise intensity V 2 ,o = 1 are simply design knobs used to determine the 
control authority. (Hence, i?o(A) = </„'(A)/?o t o a nd ^(A) = ^(AJV^o.) The system costs are 
computed assuming Vo = 0 although Vo = 0 is not assumed in the design process. This general 
philosophy is actually motivated by insights into LQG theory. However, it will suffice here to 
simply note that this philosophy was used successfully on two hardware experiments involving 
control design and implementation [Collins, Phillips, and Hyland 1991, Collins, King, Phillips and 
Hyland 1992). It should be noted that these assumptions do not influence the qualitative results 
described below. 

Below, we will compare various algorithm options. In particular, we desire to illustrate the types 
of convergence that are somet imes achieved when various bases are used to represent the controller, 
and the speed of the algorithm when various prediction options are used. We will also investigate 
what type of convergence and speed are achieved when Z/^ 1 , the inverse of the Hessian of the cost 
is not computed explicitly but is estimated using a Quasi-Newton method. The comparisons are 
all based on a MATLAB implementation of the algorithm and the program in each case was run 
on a 486, 33 MHz PC. 

We choose to base the comparison on the design of an 8th order controller (for the 8th order 
design plant). Of course, we can solve for optimal full-order controllers using Riccati equations but 
we choose this order controller because experientially we have seen that the higher the order of 
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the controller the more the algorithm struggles when a particular basis is chosen for the controller. 
Hence, we are essentially basing our comparisons on the worst-case controller order for this par- 
ticular design model. The controller that is used to initialize the algorithm is the LQG controller 
corresponding to the choice qi = 1. The algorithm is used to deform this controller into the higher 
authority controller corresponding to <72 = 0.1. 

Table 5.1 shows a comparison of the algorithm when various bases are chosen for the controller. 
Linear prediction is used in each case. In fact, it was seen experimentally that if cubic spline 
prediction was used, the algorithm performance degraded if an over-parameterized controller basis 
(i.e., tridiagonal basis or no basis) was used. This phenomenon is almost certainly due to the 
fact that in these cases the tangent vectors (0p(A)) are only estimated using (4.71) and hence are 
not accurate. As evidenced from Table 5.1, the performance of the controllable canonical form 
was worse in terms of clock time and minimum and maximum step size. The minimum step size 
of 7.8e-16 indicates substantial ill-conditioning along the homotopy path. For this example, the 
second-order polynomial form required the least number of flops although it did require slightly 
more clock time than the tridiagonaJ basis. In terms of minimum and maximum step size, the 
choice of no controller basis was better conditioned than restriction to any of the bases. 


Controller 

Basis 

Megaflops 

Real Time 
(sec.) 

No. Hessian 
Calculations 

Minimum 
Step Size 

Maximum 
Step Size 

None 

1098 

1098.2 

47 

0.01 

0.32 

Tridiagonal 

590 

880.7 

120 

0.0003 

0.08 

SPF 

51S 

930.4 

283 

0.0001 

0.04 

CCF 

828 

1524.7 

4G1 

7.8e-16 

0.02 


Table 5.1. Comparison of Controller Basis Options 


Table 5.2 shows a comparison of the algorithm when the second-order polynomial form was 
chosen for the controller and various prediction options were used. Notice that in terms of real 
time, linear prediction required only 17. S% of the time required when no prediction was used. Cubic 
spline prediction required only 5.G% of the time required when no prediction was used. The ability 
to predict along the curve described by the changing parameters is one of the practical benefits of 
formulating an optimization problem formally in terms of a homotopy. 


Prediction 

Option 

Megaflops 

Real Time 
(sec.) 

No. Hessian 
Calculations 

Minimum 
Step Size 

Maximum 
Step Size 

Correction 

Tolerance 

None 

3500 

5205.0 

1552 

6e-15 

0.01 

IHIfiHB 

Linear 

518 

930.4 

283 

1 .5e-4 

0.04 

■TSH 

Cubic 

160 

293.2 

86 

0.01 

0.08 

10“ 6 


Table 5.2. Comparison of Prediction Options for SPF Controller Basis 
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Table 5.3 shows a comparison of the algorithm when the second-order polynomial form was 
chosen for the controller, H was estimated using a Quasi-Newton (in particular BFGS) method 
and various prediction options were used. The under the Megaflop heading indicates that 
the MATLAB flop counter overflowed and so the flop data is unavailable. Notice that when the 
Quasi-Newton method was used, the prediction did not help. This is because of the inaccuracies 
in the tangent vectors due to the errors in the estimate of the inverse Hessian. Also note that by 
comparing Table 5.2 with Table 5.3, the behavior of the Quasi-Newton method was substantially 
worse than the behavior of the algorithm when the Hessian inverse was calculated exactly. In fact 
the best clock time for the Quasi-Newton method was 27 times slower than the best clock time 
when the inverse Hessian was calculated exactly. 


Prediction 

Option 

Megaflops 

Real Time 
(sec.) 

Minimum 
Step Size 

Maximum 
Step Size 

None 

* 

7960.3 

1.0e-14 

0.01 

Linear 

* 

8011.4 

1.0e-14 

0.01 

Cubic 

♦ 

8902.1 

1.0e-14 

0.01 


Table 5.3. Comparison of Prediction Options for SPF Controller Basis 
with Quasi-Newton Approximation to Inverse Hessian 


Figures 5.2 through 5.4 consider respectively the design of 2nd, 4th and 6th order controllers for 
authority levels corresponding to r/o G (1, 0.1, 0.01,. ..l.Oe — 6) and compare the optimal curves for 
an LQG controller, a reduced-order controller obtained by balancing and an optimal reduced-order 
controller. In each case, the optimal reduced-order controller performs substantially better than the 
balanced controller as the authority level is increased (i.e., is decreased). At low authority, the 
cost curves of the balanced and optimal controllers coincided, indicating that the two controllers 
are probably very similar. In fact the low authority balanced controllers were used to initialize the 
homotopy algorithm in the design of the optimal reduced order controllers as discussed in Option 
3 of Subsection 4.4. Figure 5.5 compares the optimal controllers of various orders. This type of 
figure can be used in practice to determine the order of the controller to be implemented. 

6.0 Conclusions 

This paper has presented a new homotopy algorithm for the design of H 2 optimal reduced-order 
controllers. The example of the previous section illustrated some of the features of the various al- 
gorithm options. For the test case considered, the option of estimating the inverse Hessian (7/” 1 ) 
via a Quasi-Newton method performed considerably worse than the option of actually comput- 
ing the Hessian and inverting it. The results also show the ill-conditioning that can occur when 
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a particular basis is chosen for the controller. For example, the second-order polynomial form 
was particularly ill-conditioned for the test case. In addition, the tridiagonal basis, which over- 
parameterizes the compensator, actually outperformed the second-order polynomial form in terms 
of clock time required. 

This ill-conditioning is not new. It is well known that restriction to a particular controller 
basis can cause numerical ill-conditioning or even instability (Kuhn and Schmidt 1987, Ge, Collins, 
Watson, and Davis). At least two solutions are possible. One is to have a family of minimal 
controller bases and have the algorithm switch to the basis that is best conditioned (Kuhn and 
Schmidt 1987, Ge, Collins, Watson, and Davis). Besides the second-order polynomial form and the 
controller canonical form mentioned here, another basis that could be included in this family is 
the input normal Riccati basis of (Davis, Collins, and Hodel 1992). As observed here, one can 
also use a slightly over parameterized controller basis such as the tridiagonal form. However, even 
these bases will not always be well-conditioned. One other option is to augment the cost function 
with a term that includes the squares of the free controller elements (Kuhn and Schmidt 1987). 
Unfortunately, this alternative requires a cost function that is not well motivated physically. In our 
opinion, finding practical solutions to ill-conditioning is the fundamental problem in the numerical 
computation of optimal reduced-order controllers. 


24 











control cost 


Figure 5.3. Performance Curves for the 4 th Order Controllers 
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Appendix A: Closed-Loop Matrix Derivatives 


In this appendix we present explicit expressions for the derivatives j%j, 


4t where 


e = 


vec(.A c ) 
vec (B c ) , 
vec(C c ) 


SA nr»(| 

a*’ 3110 


(A.l) 


- _ f A -BC c 
4- [fl c C A c — B c DC c \ ’ 


(A.2) 


= 



(A.3) 


where 


and 


where 


-fin — Ri 

(A.4a) 

R\i — —RiiC c 

(A.46) 

R22 = cj R?C c , 

(A.4c) 

V=\% ?»1 
, v n v 22 . 

(A.5) 

V'11 = v, 

(A. 6a) 

V12 = V u Bj 

(A. 66) 

v 22 = Bjv 2 Bj. 

(A.6c) 


It is assumed that the plant matrices (A,B,C,D), the cost weighting matrices (Ri,Ru,R 2 ) 
and the disturbance matrices (l / i , Vjj, Vj) are the following functions of A. 


A( A) B( A) 
C( A) D{ A) 


Ao 
C 0 



Bj _ Ao 
Dj Co 



(A.7) 


where 


i?i(A) tf, 3 (A) 

/? 2 ( A) 


= Lr(A)LT (A) 


(A. 8a) 


L/?(A) = L/1,0 + A (Lrj - Lfifi) 


(A.8 6) 
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and Ljifl and Lrj satisfy 


where 


•£'K,0-£'fl,0 — 


V,(A) V 12 (A) 


#1,0 

#12,0 


#12,0 

#2,0 _ 


#1,/ 

#12./ 


,/ 

#2,/ . 


] = L v (X)Ll 

(A) 


(A.8c) 

(A.8rf) 

(A.9a) 


Lv{ A) = Ly,o + A (Lvj - Lv, o) 


(A.94) 


Below, we use the notation 


Note that from (A.7)-(A.9) 


where 


where 


Mi 


dM 
d\ ' 


(A. 10) 


A tfl 
C D ~ 

A/ — Aq B j — Bo 1 
Cj — Cq Dj — Do J 

(A.ll) 

#1 #12 
_#?2 #2 

= LrLr + LrLr 

(A. 12a) 

Lr ■ 

= £/*,/ - £/?,o 

(A. 126) 

■ V 1 v 12 ' 

.v3 fcj 

j • T 

= LyL\f + LvLv 

(B.12e) 

Lv : 

= Lv,/ — Lv.o- 

(A.12d) 


The derivations of the expression for and are primarily based on the application 

of the following derivative formulas. It is assumed that A’ is an m x n matrix and A is an n X p 
matrix. 
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Derivative Formulas 


jLxA = e^A{j,:) 

(A. 14) 

-/-AX = A(:,t)e</> 
dx{j 

(A.15) 

^-,X T A = 

(A.16) 

£-AX* = Al:.i)4> 

(A.17) 

4- AXB = A{:,i)B(j,:) 

dxij 

(A. 18) 

4~AX T B = A(:J)B(i,:) 
ax a 

(A.19) 


Derivatives with respect to Bj for Oj = a Ct kt 


dA 

da c ,kt 


'0 

0 


dR 

dV 

dl> c ,kt 


= 0 
= 0 


0 




(A.20) 

(A.21) 
(A. 22) 


Derivatives with respect to Bj for 6 } = b Ci ki 


dA 

db c ,kt 

dR 

db c ,kt 

dV 

db c ,kt 




0 


0 

SYM e 


[l' l3 (:,f)-BD c K 2 (:,0)eK )T 

+ B c V,(:J)c ( n k c ] 


T 


(A.23) 

(A.24) 

(A.25) 
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Derivatives with respect to 8j for 0j = c c ^i 


where 


where 


dA 

’o -B(:,k)en} T ' 


dc c ,ki 

0 -B c D(:,k)e£ T . 


dR 

0 [-Ai2(:,*) + C T r>7i?2(:,*)]e ( „i ) 

dc C 'ki 

_SYM CjRj(:, k)e ( n} T + e ( ^R 2 {k, :)C e 

* =« 

1 



(A. 26) 
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Appendix B: Efficient Computation of H$ 

In this appendix we show how to efficiently compute Ho, using (4.26) with (4.37)-(4.39), 
(4.40)-(4.48), (4.31)-(4.33) and (A.20)-(A.28). First, we assume that $ transforms A € IR nc/Xn *«/ 
to either, complex modal form or complex Schur form, such that 

4' _1 A4' = A (B.l) 

where A € <C n * ,x " c ' is diagonal or upper triangular. The pre- and post-multiplying (4.31) respec- 
tively by and '4', pre- and post-multiplying (4.32) by respectively 4> -1 and 9~ H , and pre- 
post-multiplying (4.33) by T~ l and T give 

0 = A *P (i) + P U) A + (A (i)T P + PA U) + R {j) ) (B.2) 

0 = AQ U) + Q U) A* -I- (A U) Q + QA u)T + V ( ») (B.3) 

Z U) - QU)p + QpU) (B.4) 

where 

pU) — pU)\Jf 
qU) _ y-iQU)y-H 

Z U) = 4 » ~'Z U) V 
P= 4- H PV 
Q = 4 '- 1 04 » -// 

A u) = 9~ l A U) 9 

pU) _ 

\Ai) — 

Next, partition 4^ as 

4' = , 4», € IR n * x " £ ', 4»2 € IR n ‘ Xflc< (B.13) 

and partition 4' -1 as 

4' -1 = [(4'- 1 )i,( 4 ,-1 )2], (4<- , ) 1 e IR n ‘ ,>(r, ^ (r'jjer-'*"'. (B. 14 ) 


(B.5) 
(B.6) 
(B.7) 
(B.8) 
(B.9) 
(B.10) 
(B-11) 
(B. 12) 
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Also, define 


b = (y-')iB 

(B.15) 

C = C9 1 

(B.16) 

Rn = 9?Rn 

(B.17) 

Vl2 = (»- 1 ) 1 V 12 

(B.18) 

B e = (9- l hBc 

(B.19) 

C c = C c * 2 . 

(B.20) 


Now, recall that where 0j represents either a c .fc*,&c,JM,c c >*- It then follows from 

(B.10)-{B.12) and (A.20)-(A.28) that A^\ R^ and are given as follows. 


for 0j = a c ,k< 

= :) (B.21) 

= 0 (B.22) 

V U) = 0 (B.23) 

for 8j = b c , kf 

i 0) = *; 1 {:,k)[C{t,:) + D(l,:)C e ] (B.24) 

= 0 (B.25) 

V {}) = {[V r 12 (:,£)-BZ) c V 2 (:,0 + 5 c V 2 (:,/)][(»- 1 ) 2 (:,A-)] H } 

+ {[V 12 (:,f) + B c V 2 (:,0][(^ _1 )2(:,A:)] f/ } . (B.26) 

for 8j = c c , k / 

A 0) = [£(:,*) + £ C I?(:,*)]¥ 2 (£,:) (B.27) 

R U) = { [(* _1 ) 2 (/, 0] H [~Rn(:,k) H + R 2 (k,:)D c C + R 2 (k, :)C c ] } 

+ { [(*' ■' M* , : )] H [-*»(:, k) H + R 2 (k, :)C e ] } (B.28) 

V U) = 0 (B.29) 
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Note that (B.21)-(B.29) allow the transformation of R W and to the modal 

basis to be performance very efficiently. 

Now, it follows from (B.5)-(B.7) that 

pU) — ty-H 

QU) - 

Z U) = 


or, equivalently, 


p(i) 

Ml 

p(i) ■ 
M2 


'($- 1 ){ f pW(^' 1 ), (»' 


6(i) T 

M2 

p(» 

^22 J 


(*■ 

-I)ffp0)($-X) 2 


$<,?• 





-Q[f 

«?. 





■«» 

7( j) " 




*U) 

L ^21 

v(j) 
*22 J 


’M<»(<p - , ) 1 <M (j) (*- 1 ) 2 



It follows from (4.37), (B.35), and (B.4) that 

H‘ Ac (Z U) ) = H ( i] = ViQMpW-'h + 2 

It follows from (4.38), (B.33), (B.35), (B.13), (B.14), and (B.16)-(B.18) that 

± H ( i] = 2 [*?(P {j) M HB c) + P2(Q {j) C HB c)] 

where 

A = * 2 H p 

Mhbc = A 2 + B c Vi + QCj/Bc 
Chbc = C - DC c . 

Similarly, it follows from (4.39), (B.34), (B.35), (B.13)-(B.15), (B.17) and (B.1S) that 

= H% = 2 [(M H ccQ {3) )*f - ( B%ccP U) )Q 2 ] 


of Schur 


(fl.30) 

( 5 . 31 ) 

(B.32) 

(B.33) 

(B.34) 

(B.35) 

(B.36) 

(B.37) 


(B.38) 

(B.39) 

(B.40) 


(B.41) 
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where 


Ql = * 2 Q 

Mhcc = RiC c - - BhccP 

Bhcc = B + B C D. 


Finally, substituting (B.33)-(B.35) into (4.40)-(4.48), using (B.13)-(B.18) and 
definitions (B.36), (B.37) and (B.41) gives the following. 

Derivations with respect to a Ct ki 


Ma. _ ff U) 

da cM ~ 


QHbc 

d&cM 

dHc. 

d&cM 




Derivations with respect to b Ct kt 


_ Mi) 

db cM A ‘ 

= H b1 - 2(A(» _, )2 (:,*)) W,:) 

OVc,kt 

= H l c *} - 2D(f.,:) T {P 7 (k,:)Q?) 

Obc.kt 

Derivations with respect to c Ci kt 


dH^ 

@Cc,kt 



BHb, 

dcc.ke 

dHc, 

dc c ,k( 




lihW, :)")£>(:, *) T 




(B.42) 

(B.43) 

(B.44) 

recalling the 

(£. 45 ) 

(5.46) 

(5.47) 

(5.48) 
(B.49) 

(5.50) 

(5.51) 
(B.52) 
(5.53) 
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Abstract 

It has been observed numerically that low authority LQG compensators are often nearly non- 
minimal. However, to date a rigorous justification for this phenomenon has not yet been established. 
This paper helps provide the needed theoretical foundation. In particular, it is shown that for both 
continuous-time and discrete-time stable systems, by proper choice of the structure of the design 
weights, the corresponding LQG compensator becomes nonminimal as the control authority is de- 
creased. Thus, the results provide a partial explanation of why the suboptimal controller reduction 
methods tend to work best at low control authority. The results also can be used as rigorous guide- 
lines to efficiently initialize homotopy algorithms for directly synthesizing optimal reduced-order 
controllers. The restriction to stable systems is not necessarily limiting since the freedom involved 
in defining a homotopy allows this assumption to always be satisfied. 
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1. Introduction 


The development of linear-quadratic-gaussian (LQG) theory [1*3] was a major breakthrough in 
modern control theory since it provides a systematic way to synthesize high performance controllers 
for nominal models of complex, multi-input multi-output systems. However, one of the well known 
deficiencies of an LQG compensator is that its minimal dimension is usually equal to the dimension 
of the design plant. This has led to the development of techniques to directly synthesize optimal 
reduced-order controllers [4—17] and techniques to synthesize reduced-order approximations of the 
optimal full-order compensator (i.e., controller reduction methods) [18-23]. 

The controller reduction methods almost always yield suboptimal (and sometimes destabilizing) 
reduced-order control laws since an optimal reduced-order controller is not usually a direct function 
of the parameters used to compute or describe the optimal full-order controller. Nevertheless, these 
methods are computationally inexpensive and sometimes do yield high performing and even nearly 
optimal control laws. An observation that holds true about most of these methods is that they 
tend to work best at low control authority [17, 21, 23]. However, to date no rigorous explanation 
has been presented to explain this phenomenon. 

One of the purposes of this paper is to provide a partial explanation as to why the suboptimal 
projection methods tend to work at low control authority. The discussion here focuses on stable 
systems. It is shown that if the state weighting matrix R\ or disturbance intensity (or covariance 
for discrete systems) Vj has a specific structure in a basis in which the A matrix is upper or lower 
block triangular, respectively, then at low control authority the corresponding LQG compensator 
is nearly nonminimal and can hence be easily reduced to a nearly optimal reduced-order controller. 
The conditions presented for R\ and V\ often are satisfied or nearly satisfied in practice. Hence, for 
stable systems the results proved in this paper do offer one explanation of why suboptimal controller 
reduction methods often provide nearly optimal control laws at low authority. The results can also 
be used as guidelines for choosing R\ and \\ such that suboptimal controller reduction methods 
yield “good 1 ' reduced-order controllers. 

Suboptimal controller reduction methods can be used to initialize algorithms for synthesizing 
optimal reduced-order controllers. Of particular interest are the homotopy algorithms of [11. 15-17] 
since they are based on allowing the plant and weights defining an optimization problem to vary 
as a function of the homotopy parameter A € [0.1]. These homotopy algorithms rely on choosing 
the initial plant and weights so that the corresponding LQG compensator is easily reduced to a 
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nearly optimal reduced-order compensator of the desired dimension. Hence, the results presented 
here provide some rigorous guidelines for initializing these algorithms. Note that the restriction to 
8table systems is not necessarily limiting since the freedom involved in defining a homotopy allows 
this assumption to be satisfied. However, future work will focus on theory that directly applies to 
unstable systems. 


Notation 

IR,IR rXi ,IR r 

IE 

X > 0,X > 0 

Orxn^r 

Ir 

vec(-) 


real numbers, rxs real matrices, IR rxl 
expected value 

matrix X is nonnegative definite, X is positive definite 
rxs zero matrix, r x r zero matrix 
r x r identity matrix 

the invertible linear operator defined such that 

vec 5 = [«7 ‘ ' ■ ,S J] T > § € 

where Sj € IR P denotes the j th column of 5. 


2. Low Authority LQG Compensation: Continuous-Time Systems 

Consider the n t/l -order Unear time-invariant plant 

i(t) = Ax(t) + Bu(t) + Diw(t), (2.1a) 

y(t) = Cx(t) + D 2 w(t), (2.1 b) 

where (A, B) is stabiUzable, {A,C) is detectable, x € IR n ,u € IR m ,y € IR/, and w € IR rf is a 
standard white noise disturbance with intensity Id and rank D 2 — l. The intensities of D\w(t) and 
D 2 w(t) are thus given, respectively, by Vj = D\Dj > 0, and V 2 = D 2 Dj > 0. For convenience, we 
assume that V 12 = D\Dj = 0, i.e., the plant disturbance and measurement noise are uncorrelated. 
Then, the LQG compensator 


x c (t) = A c x c (0 + B c y(t), 
u(t) = —C c x c {t), 

for the plant (2.1) minimizing the steady-state quadratic performance criterion 

J(A c ,B c .C c )= lim -E f [t t {s)R\x(s) + v r (s)R'>u(s)]d$. 
i-* 1 J 

o 


(2.2a) 

(2.26) 


(2.3) 
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where Ry > 0 and R 2 > 0 are the weighting matrices for the controlled states and controller input, 
respectively, is given by: 

A c = A — YP — QY. (2.4a) 

B c = QC T V 2 7 1 , C c = R; l B r P, (2.46, c) 

where E — BR^ *5^, E = *C, and P and Q are the unique, nonnegative-definite solutions 

respectively of 

0 = A r P + PA + Ry- PY.P, (2.5) 

0 = AQ + QA t + V, - QtQ. (2.6) 

Furthermore, the “shifted” observability and controllability grammians [18, 24] of the compensator, 
P and Q, are the unique, nonnegative-definite solutions respectively of 

0 = (A - Q£) T P + P(A -Qt) + PEP, (2.7) 

0 = (A - ZP)Q + Q(A - ZPf + QXQ. (2.8) 

Although a cross-weighting term of the form 2i T (t)f? 12 u(t) can also be included in (2.3), we shall 
not do so here to facilitate the presentation. The magnitudes of R 2 and V 2 relative to the state 
weighting matrix Ry and plant disturbance intensity Vi govern the regulator and estimator au- 
thorities, respectively. The selection of R 2 and V 2 such that ||# 2 || >> ||#i||, or ||V 2 || >> flVjjJ, 
yields a low authority compensator. It has been observed numerically that low authority LQG 
compensators are often nearly nonminimal [17, 21]. This section provides a rigorous justification 
for this observation when the open-loop plant is stable and (A, Ry) or (A,Vi) have a particular 
structure. In order to prove this result, we first exploit some interesting structural properties of 
the solutions of the Riccati equations and Lyapunov equations assuming the coefficient matrix A 
and the constant driving term Ry have certain partitioned forms. 


Lemma 2.1. Suppose 
A = 


Ay 0 
A 2 i A 2 


, B = 


By 

b 2 


, Ry = 


#1,1 o 

[ 0 0„_ nr J 


(2.9a, 6,c) 


where Ay,Ry,y € IR"'*"', By € IR" rXm , R 1A > 0. 


(i) If (A, B) and (Ay, By) are stabilizable, then the unique, nonnegative-definite solution of 
the Riccati equation: 

0 = A 1 P + PA + Ry - PBB J P. 


(2.10) 



is given by 

p _ [ Pi 0 

l 0 0 B _ nr 

where the n r x n r matrix Pi is the unique, positive-definite solution of 

0 = Aj Pi + PiAi + /?i,j — P\BiBj P\. 


( 2 . 11 ) 


( 2 . 12 ) 


(»*) If A is asymptotically stable, then the unique, nonnegative-definite solution of the Lyapunov 
equation: 

0 = A t P + PA + R u (2.13) 


is given by 


p=\ Pl 0 1 

lo 0„_ n J’ 


where the n r x n r matrix P\ is the unique, positive-definite solution of 

0 = A{Pi + P\Ai + 


(2.14) 


(2.15) 


Proof. 


(i) Since (A, B) is stabilizable and R\ > 0, it follows from Theorem 12.2 of [25] that there 
exists a unique, nonnegative-definite solution of the Riccati equation (2.10). Similarly, the 
assumptions that (j4i,jBi) is stabilizable and > 0 imply that there exists a positive- 
definite matrix P\ satisfying the Riccati equation (2.12). Using (2.12), it follows by con- 
struction that (2.11) is the solution of (2.10). 

(it) This is a special case of the Riccati equation of property (i). □ 


The following lemma states the dual of Lemma 2.1 if the coefficient matrix A is upper block 
triangular and V\ is upper block diagonal. 


Lemma 2.2. Suppose 


A = 


A\ A\2 
0 A 2 


, c = [Cl C 2 ], Vi = 


Vi,i 0 

o 0 „_„, 


(2.16a, 6,c) 


where A\,Vi^ € IR" rX " r , V'i.i > 0. 


(i) If {A,C) and (A\,C\) are detectable, then the unique, nonnegative-definite solution of the 
Riccati equation: 

0 = AQ + QA 1 + \\ - QC t CQ, ( 2 . 17 ) 
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is given by 



0 

On — n r 


y 


where the n r x n r matrix Q\ is the unique, positive-definite solution of 


0 = AiQi + QiAj + Vi,, - QiCj CiQi. 


(2.18) 


(2.19) 


(n) If A is asymptotically stable, then the unique, nonnegative-definite solution of the Lyapunov 
equation: 

0 = AQ + QA t + Vi, (2.20) 


is given by 


Q = 


Q i o 

|_ 0 0 n _ n , 


( 2 . 21 ) 


where the n r x n r matrix Q\ is the unique, positive-definite solution of 


0 = A\Q\ + Q\A l + V,,,. 


( 2 . 22 ) 


Proof. The proof is dual to the proof of Lemma 2.1. □ 

The following theorem shows that, with proper choice of the weighting matrices, a low authority 
LQG controller for a stable plant is nearly nonminimal. The proof of this theorem relies on the 
above two lemmas. 


Theorem 2.1. Consider the plant given by (2.1). 
(i) Suppose 


A = 


A x 0 

An Ai 


, Ri = 


Ri'i 0 

0 0 n _ n , J ’ 


(2.23a, 6) 


where j4i , ^2, ,x € nt nrXnr , i? ltl > 0, and A is asymptotically stable. Let 


V 2 = f3V 2 


(2.24) 


where V 2 is some finite, positive-definite matrix and /? € IR is a positive scalar. Then 

lim rank (QP) < lim rank P < n r , (2.25) 

P—+OC — OG 

where Q and P are the shifted controllability and observability grammians of the cor- 
responding LQG compensator, satisfying (2.8) and (2.7). respectively. Equivalently, for 



6 > 0, there exists N such that for all 0 > N, A„ r+1 < tfA nr , where A * represents the i th 
eigenvalue of QP and Aj > A 2 > ... > Aj > Aj + i... > 0. 


( ii) Suppose 


_U a 12 ] . fv,., 0 1 

[o A 2 J ’ 1 ~ l 0 o n _„ r J ’ 


(2.26a, 6) 


where A\, Vi 4 € IR" r " r , Vi 4 > 0, and A is asymptotically stable. Let 


R2 = aRi, 


(2.27) 


where iZ 2 is some finite, positive-definite matrix and a € IR is a positive scalar. Then 


lim rank ( QP ) < lim rank Q < n r , (2.28) 

Or —*■ 00 a— ♦ 00 

where Q and P are the shifted controllability and observability grammians of the cor- 
responding LQG compensator, satisfying (2.8) and (2.7), respectively. Equivalently, for 
6 > 0, there exists N such that for all a > N, A nr +j < 6A„ r , where A; represents the i th 
eigenvalue of QP and A! > A 2 > ... > A, > A,- + i... > 0. 


Proof. 

* B " 2 S 1 

(i) Partition B = ^ and E = ^ ^ , conformal to A in (2.23). The assumptions 

(2.23) and that A is asymptotically stable imply that (.4, B ) and (4 1? B\) are both stabiliz- 
able. Thus, it follows from property (t) of Lemma 2.1 that the unique, nonnegative-definite 
solution P of the Riccati equation (2.5) has the structure given by (2.11), which implies 
that 

O' 


PEP = 


fiEiPi 

0 


0 


(2.29) 


Thus, noting the special partitioned structures in (2.29) and (2.23), and that A is asymp- 
totically stable, it follows from property (ii) of Lemma 2.1 that there exists 


Po = 


Pi 

0 


0 

n r 


which is the unique, nonnegative-definite solution of 


0 = A r P 0 + P 0 A + PEP , 


(2.30) 


(2.31) 


where n T x n r matrix P\ is the unique, nonnegative-definite solution of 


Next, computing (2.31) - (2.7) and using (2.24), yields the following modified Lyapunov 
equation: 

0 = A T AP + A PA + /?- 1 [(C t V 2 " 1 CQP) + (C T Vf'CQP) T ). (2.32) 


where 

AP = P 0 - P. 


(2.33) 


Since A is asymptotically stable and Q and P satisfy (2.6) and (2.7), respectively, Q and 
P are bounded for all /3. Thus, (2.32) implies that lim/j_oo ||A.P|| = 0. Hence, for c > 0, 
there exists M such that for all /? > M, ||A.P|| < e. Using (2.33), it follows that 


lim P = Po 

0-+OO 


A o" 

0 0 • 


(2.34) 


Thus, limbec rank ( QP ) < lim^—oo rank P = n r , which implies the following inequaE- 
ties of the eigenvalues of QP. Suppose A t represents an eigenvalue of QP and Ai > A 2 ... > 
A t > Aj + i... > 0. Then, for 6 > 0, there exists jV such that for all (i > A r , A nr + 1 < 6X nr . 


(«) The proof is dual to the proof of (i). □ 

Remark 2.1. Theorem 2.1 provides two ways of weighting matrices selection resulting in a 
nearly nonminimal, low authority LQG compensator for a stable plant. The first approach starts 
by transforming the plant A into coordinates such that A has the representation as in equation 
(2.23a) after transformation. Then select the weighting matrix J?i with the partitioned form as in 
(2.23b) and with rank R\ = n r . By decreasing the authority of the compensator, or. equivalently r 
increasing 1 1 1 ^ 1 1 or /?, the LQG compensator approaches nonminimality with minimal dimension of 
n r . Using a dual approach, with A and Vi partitioned as in (2.26), by increasing ||/?2 II or a , 
resulting LQG compensator approaches nonminimality. 


Remark 2.2. Note that if A is in a modal form, then it satisfies both (2.23a) and (2.26a) of 
Theorem 2.1. In this case, R\ given by (2.23b), describes a state weighting matrix in which only 
the states pertaining to selected modes are weighted. Similarly, V\ given by (2.26b) describes a 
disturbance that excites only certain modes. It is not uncommon for these conditions to be satisfied 
or nearly satisfied in practice. 


Remark 2.3. The suboptimal controller reduction methods of [18-23] characterize the redu- 
ced-order controller by a projection or some other type of reduction of the LQG controller. It 
has been observed that these suboptimal reduced-order controllers for the low-authority control 



problem will yield virtually the same cost as the LQG controller. According to Theorem 2.1, for 
a stable plant and with proper choice of the weighting matrices, the LQG controller for a low 
authority control problem is nearly nonminimal, which provides a theoretical justification for the 
above observation. 


Remark 2.4. The homotopy algorithms for reduced-order dynamic compensation problems 
developed in [15-17] are based on allowing the plant and weights defining an optimization problem 
to vary as functions of the homotopy parameter A € [0, 1]. In particular, it is assumed that 


A( A) B{ A) 
C(A) 0 


Ao Bq 
C o 0 


+ K 


[% « ] - [ 


Aq Bo 

Co 0 


). 


where 

and and Z/*j satisfy 

£k,oLh,o = 

and 

where 

and Xv'.o and Lvj satisfy 

LvfiLlo = 


■fti(A) R i2 {X) _ Tf\\rT(\\ 
Rj 2 ( A) R 2 ( A) J " L *W L i&Ah 


Lr(X) = Lfifl + A (Lrj - Lr^o), 


R i,o 
■ft 12,0 


R 12,0 

^2,0 


L R,f L R,J - \^T ' f 


Vi(A) V n ( A) 
IVS(A) V 2 (X)\ 


= L v {X)Ll(X), 


R 12 ,f 

R 2J J 


Lv{ A) = Lv,o -I- A (Lv,j — fv,o)i 


V, f0 

t'12,0 


Vi2fi 
V 2 ,o j 


LvjLyj — 


ViJ 

^12./ 


Via ,/ 

VtJ 


Note that the above equations imply that .4(0) = Ao, B{ 0) = Bo, etc ... which are the ini- 
tial set of system matrices and that 4(1) =4/, B( 1) = Bj, etc ... which are the final and 
given system matrices, To initialize the homotopy algorithm efficiently, the designer can choose 
( Ao, Bo, Co, Rifi, R\ 2 ,o, r 2 ,o, Vj,o> k’i 2 ,o • 1'2 .o ) • to correspond to a low authority control problem with 
stable open-loop plant as stated in Theorem 2.1, for which a nearly optimal reduced-order controller 
may be easily obtained by balanced controller reduction [18] or an alternative suboptimal controller 
reduction method [19-23]. 
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3. Low Authority LQG Compensation: Discrete-Time Systems 


In this section, we consider the discrete- time counterpart of the previous section. In particular, 
a rigorous justification is provided for a nearly nonminimal low authority discrete-time LQG com- 
pensator when the open-loop plant is stable and certain weighting matrix has specific structure. 
Consider the n </l -order linear, discrete time-invariant plant 

x(k + 1) = Ax{k ) + Bu(k) + D x w(k), (3.1a) 

y(k) = Cx(k) + D 2 w(k), (3.1b) 

where (A,B) is stabilizable, ( A,C ) is detectable, x € IR n ,u € lR m ,y € and w € IR d is a 
standard white noise disturbance with covariance Id and rank D 2 = /. The covariances of D\w(k) 
and D 2 w(k) are thus given, respectively, by Vi = D x Dj > 0, and V 2 = D 2 D 2 > 0- For convenience, 
once again we assume that Vj 2 = D\Dj = 0. Then, the LQG compensator 

x c (k + 1) = A c x c (k) + B c y(k), (3.2a) 

u(k) = -C c x c (k) - D c y(k), (3.2 b) 

for the plant (3.1) minimizing the steady-state quadratic performance criterion 

J{A c ,B c ,C c ,D c ) = lim IE[x T (ib)i? 1 a:(Jt)-|-u T (it)ii 2 u(l:)], (3.3) 

fc— • OO 

where R\ > 0 and Ri > 0 are respectively the weighting matrices for the controlled states and 
controller input, is given by [26]: 


A c — A — Q a^ 2 a C “ BR 2 a Pdf (3.4a) 

B c = Q a V 2 -\ C c =R^P a . Z) c = R 2a 1 B T PAQC T V 2 -\ (3.46, c,d) 

where 

Q a = AQC J , P a = B 1 PA. V 2a = V 2 + CQC r , R 2a = R 2 + B T PB, (3.5. 6, 7.8) 

and P and Q are the unique, nonnegative-definite solutions respectively of 

P = A 1 PA + 7?, - Pj R 2 g P a , ( 3.9 ) 

Q = AQA 1 + V, - QaVfa'Ql- (3.10) 

<> 



Furthermore, the “shifted” observability and controllability grammians of the compensator, P and 
Q, satisfy 


P = {A- QaV{ a l C) J P{A - QaV^C) + (P a - R 2 a D c C) T (P a - R 2a D c C ). (3.11) 

Q = (A- BR£Pa)Q(A - BR^P a f + (Q a - BD c V u )V^{Q t - BD e V 2a f , (3.12) 

and P and Q are nonnegative definite. 

As in the continuous-time case a cross- weighting term of the form 2x T (k)R\ 2 u(k) can also be in- 
cluded in (3.3), we shall not do so here to facilitate the presentation. Similar to the continuous-time 
compensation problem, the magnitudes of R 2 and V 2 relatively to iZj and V\ govern the regulator 
and estimator authority, respectively. The following theorem is the discrete-time counterpart of 
the continuous-time result stated in Theorem 2.1. It provides a rigorous justification for a nearly 
nonminimal low authority discrete-time LQG compensator when the open-loop plant is stable and 
Ri or \\ has certain structure. 


Theorem 3.1. Consider the plant described in (3.1). 
(i) Suppose 


A = 


A, 0 

A 2 i A 2 


, /?i = 


Ri 

0 


.1 0 ] 

On-nJ’ 


(3.13a, 6) 


where A\,R\ t \ € IR nrXrir , f?i t i > 0, and A is asymptotically stable. Let 


V 2 = pv 2 . 


(3.14) 


where V 2 is some finite, positive-definite matrix and /? € IR is a positive scalar. Then 


lim rank (QP) < lim rank P < n r , (3.15) 

/3— *>oo (3—>oo 

where Q and P are the shifted controllability and observability grammians of the corre- 
sponding LQG compensator, satisfying (3.11) and (3.12), respectively. Equivalently, for 
b > 0. there exists N such that for all /? > N, A„ r+ i < £A„ r , where A, represents the i th 
eigenvalue of QP and Ai > At > ... > X, > A l+ i... > 0. 

(it) Suppose 

. _ A] -4]2 
“ 0 A 2 


li = 


Vt.i 

0 


0 

On — n. 


i3.16«. 6) 
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where .4i. V\ t \ € nt' 1 '*' 1 '', V x l > 0, and A is asymptotically stable. Let 

Ri = ^R-2 i (3.17) 

where R> is some finite, positive-definite matrix and a € IR is a positive scalar. Then 

lira rank ( QP ) < lim rank Q < n r , (3.18) 

a— *oo a— *oo v ' 

where Q and P are the shifted observability and controllability grammians of the corre- 
sponding LQG compensator, satisfying (3.11) and (3.12), respectively. Equivalently, for 
b > 0, there exists N such that for all a > JV, A nr+ i < SX nr , where A,- represents the i th 
eigenvalue of QP and A t > A 2 > ... > Aj > Aj +1 ... > 0. 

Proof. The proof is similar to the proof of Theorem 2.1. □ 


4. Numerical Illustrative Examples 

To illustrate the proper choices of the weighting matrices resulting in a nearly nonminimal, 
low authority LQG compensator for a stable continuous- time plant, consider a simply supported 
beam with two collocated sensor /actuator pairs. Assuming the beam has length 2 and that the 
sensor/actuator pairs are placed at coordinates a = , and b = ||, a continuous-time model 

retaining the first five modes is obtained: 

x = Ax + Bu + D\w, y = Cx + D 210 , 


where 


f 0 

1 

[ 

1 1 

0 

1 1 

[ 

0 1 1 

l-l - 

0.01 

’ [-16 

- 0 . 04 J ’ 

’ [-81 

- 0 . 09 J 1 

l- 

■256 — 0.16 J 

- 0.2174 

0 

0.4245 

0 - 0.6112 0 

0.7686 

0 

— 0 . 8893 "! T 

- 0.8439 

0 

- 0.9054 

0 - 0.1275 0 

0.7686 

0 

0.9522 J ’ 


The noise intensities are Vj = D\Dj = O.l/io and V 2 = D 2 DJ = Ph, and it is assumed that Vu = 
D\Dj = 0. The design objective is to minimize the continuous- time cost J = lim ( _oo IE[x T i2ix + 
tt T i? 2 u], where R 2 — al 2 . Note that the magnitude of the positive real numbers a and /? are the 
indicators of the controller authority level. For this particular plant, A has the representation as 
in (2.23a) and (2.26a) with Au = 0 and A 21 = 0 , respectively. Here, we illustrate the results of 
property ( t ) of Theorem 2.1 for the cases of n r = 2 and n r = 6. Setting a = 0.1, by selecting the 
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weighting matrix R\ = ^ jj , and increasing 0 (hence, decreasing the compensator authority), 

the resulting LQG compensator approaches nonminimality with minimal dimension of n r or, equiv- 
alently, X ' r + ‘jS > ^ ) — ► 0 where A, is the sorted (in descending order) i th eigenvalue of QP. Figure 1 
shows the ratio curve for n r = 2 with 0 G (0.01,0.1, 1, 10, 10 2 , 10 3 , 10 4 , 10 s , 10 5 6 ). The curve clearly 
indicates that the ratio decreases as 0 increases. To illustrate that suboptimal controller reduction 
methods yield nearly optimal reduced-order compensators for low authority control problems, Fig- 
ure 1 also shows the norm of the cost gradient of the 2 n<i -order controller obtained by balancing. 
The cost gradient is defined as [(vec (vec -§^) T (vec ] T - The cost gradient curve 

indicates the balanced controller approaches the optimal reduced-order compensator as 0 increases, 
or as the control authority decreases. Figure 2 shows the eigenvalue ratio of the LQG controller for 
n r = 6 and the norm of the cost gradient of the corresponding 6 tfc -order balanced controller. 

Conversely, if the weighting term for the above example does not have the structure given by 
(2.23b), decreasing the controller authority (i.e., increasing 0) may not yield a nearly nonminimal 
LQG compensator. As a result, the norm of the cost gradient of the corresponding 2 nd -order 
balanced controller does not approach zero as the control authority decreases. This is illustrated 
in Figure 3 for n r = 2 and R\ = I\ o- Note that for this particular example, at 0 = 0.01 the 
balanced controller destablizes the closed-loop system and hence the norm of the cost gradient 
becomes infinite. 

5. Conclusion 

By exploiting structural properties of the solutions of the Riccati equations and Lyapunov 
equations, this paper shows that for both continuous-time and discrete-time stable systems, if the 
coefficient matrix A and driving weighting term R\ (or Vj) have specific structures, the corre- 
sponding LQG compensator becomes nonminimal as the control authority is decreased. This result 
provides a partial explanation of why suboptimal projection methods tend to work best at low 
authority. This paper also establishes some rigorous guidelines to initialize homotopy algorithms 
for directly synthesizing optimal reduced-order controllers. In particular, to initialize the homotopy 
algorithm efficiently the designer can choose the plant and weighting matrices to correspond to a 
low authority control problem with stable open-loop systems as stated in Theorem 2.1 or 3.1. In 
this case, a nearly optimal reduced-order controller may be easily obtained using an appropriate 
suboptimal controller reduction method such as balancing since the resulting LQG controller is 
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nearly non-minimal. These results are clearly illustrated by numerical examples. 

Conversely, if the structure of the plant and weighting matrices do not satisfy the conditions 
specified in Theorem 2.1 or 3.1, the resulting LQG compensator is not necessarily nearly minimal 
even at low control authority. In this case, reduced-order controllers obtained by suboptimai 
projection methods may not be nearly optimal even at low authority. This result is illustrated in 
the last example with a reduced-order controller obtained by balancing. 
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cost gradient/eigenvali 


6 , 



Figure 2. Non- minimality indicator ( ) °f LQG controller ( — ) and the norm of 

the cost gradient of the 6 t/l -order balanced controller ( ) vs control authority (0) for n r = 6 

and “good structured” R\ 



gradient/eigenvalue ratio 



Figure 3. Non-minimality indicator of the LQG controller (— ) and the norm of 

the cost gradient of the 2 n< *-order balanced controller (- - -) vs control authority (/?) for n r = 2 
and w bad structured” iZj 
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